Instantly share code, notes, and snippets.

# richarddmorey/delta_t.R

Last active Feb 27, 2017
posterior for delta given only t and N
 ### Data t = 2 N = 20 rscale = sqrt(2)/2 ### Begin utility functions posterior = Vectorize(function(delta, t, N1, N2 = NULL, rscale = sqrt(2)/2, log = FALSE){ effN = ifelse(is.null(N2), N1, N1*N2/(N1+N2)) df = ifelse(is.null(N2), N1 - 1, N1 + N2 - 2) log.post = dt(t, df, delta * sqrt(effN), log = TRUE) + dcauchy(delta, scale = rscale, log = TRUE) ifelse(log, log.post, exp(log.post)) }, "delta") pdf.posterior = function(delta, ...){ norm.const = integrate(posterior, -Inf, Inf, ..., log = FALSE)[[1]] posterior(delta, ..., log = FALSE) / norm.const } cdf.posterior = function(delta, ...){ norm.const = integrate(posterior, -Inf, Inf, ..., log = FALSE)[[1]] x = integrate(posterior, -Inf, delta, ..., log = FALSE)[[1]] x / norm.const } invcdf.posterior = Vectorize(function(p, t, N1, N2 = NULL, rscale = sqrt(2)/2){ norm.const = integrate(posterior, -Inf, Inf, t = t, N1 = N1, N2 = N2, rscale = rscale, log = FALSE)[[1]] effN = ifelse(is.null(N2), N1, N1*N2/(N1+N2)) d.est = t / sqrt(effN) d.se = 1/sqrt(effN) limits = qnorm( pnorm( qnorm(p, d.est, d.se ) + c(-1,1)*d.se, d.est, d.se), d.est, d.se) optimize( function(delta)(cdf.posterior(delta, t = t, N1 = N1, N2 = N2, rscale = rscale) - p)^2, limits)\$minimum }, "p") ### End utility functions ### Compute interesting quantities # 95% credible interval ci = invcdf.posterior(c(.025, .975), t, N, rscale = rscale) # posterior median post.med = invcdf.posterior(.5, t, N, rscale = rscale) # posterior mean exp.delta = integrate( function(x,...) x * pdf.posterior(x, ...), -Inf, Inf, t=t, N1=N, rscale = rscale )[[1]] # posterior variance/sd exp.delta2 = integrate( function(x,...) x^2 * pdf.posterior(x, ...), -Inf, Inf, t=t, N1=N, rscale = rscale )[[1]] post.var = exp.delta2 - exp.delta^2 post.sd = sqrt(post.var) # Create plots dd = seq(-.5, 1.5, len = 100) plot(dd, pdf.posterior(dd, t, N, rscale = rscale ), ty='l', lwd=2, ylab = "Density", xlab = expression(delta), axes = FALSE) lines(dd, dcauchy(dd, scale = rscale), lty=3, col="gray") axis(1) box() abline(h=0, col="gray") abline(v = ci, col="red", lty=2) abline(v = post.med, col="darkgreen", lty=2) abline(v = exp.delta, col="blue") segments(exp.delta - post.sd, par()\$usr[4]/4 ,exp.delta + post.sd, par()\$usr[4]/4, col="blue") ## If you want it, here's a Bayes factor against delta = 0 B10 = dcauchy(0, scale = rscale) / pdf.posterior(0, t, N, rscale = rscale) B10 # test against BayesFactor library(BayesFactor) ttest.tstat(t, N, rscale = rscale, simple = TRUE)
Owner Author

### richarddmorey commented Feb 17, 2017

 @rezae, I just made a small change. Try it now.