Skip to content

Instantly share code, notes, and snippets.

@richarddmorey richarddmorey/delta_t.R
Last active Feb 27, 2017

Embed
What would you like to do?
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)
@richarddmorey

This comment has been minimized.

Copy link
Owner Author

commented Feb 17, 2017

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.