scale-priors.Rmdlibrary(rstap)
library(dplyr)
library(tidyr)
library(ggplot2)
library(purrr)
library(plot3D)Starting with a simple model
\[ E[Y] = \alpha + Z_1\delta_{sex} + X(\theta)_{FF}\beta_{FF} \] We’ll calculate the likelihood for different values of \(\theta\) and \(\beta_FF\) assuming that \(\alpha and \delta\) are given
data("homog_subject_data")
data("homog_distance_data")
homog_subject_data
homog_distance_data
alpha <- 22
delta <- -.8
sigma <- 2.3
thetas = seq(from = 0.01, to = 10, by =0.1)
betas = seq(from = 0.01, to = 10, by = 0.1)Z <- model.matrix(y ~ sex, data = homog_subject_data)
Xs <- homog_distance_data %>% crossing(thetas) %>% group_by(subj_id,thetas) %>%
summarise(exposure = sum(pracma::erfc(Distance/thetas))) %>% arrange(thetas) %>%
split(.$thetas)
Xs <- map(Xs,function(x) x %>% mutate(X_tilde = (exposure - mean(x$exposure))/sd(x$exposure)))
lls <- map(betas,function(a) map_dbl(Xs,function(x) sum(dnorm(x = homog_subject_data$y,mean = Z %*% c(alpha,delta) + x$X_tilde * a,sd = sigma,log = T))))
llmat <- Reduce(f = rbind,x = lls)
rownames(llmat) <- as.character(betas)
colnames(llmat) <- as.character(thetas)par(mfrow=c(1,2))
persp3D(x = betas, y = thetas, z = llmat,phi = .6,theta = 0.5,xlab="beta",ylab="theta",zlab="log-lik",ticktype ='detailed')
persp3D(x = betas, y = thetas, z = llmat, phi = 3, theta = 90, xlab="beta",ylab="theta",zlab="log-lik",ticktype = "detailed")data_frame(beta = betas[1],
thetas = thetas,
ll = llmat[1,]) %>%
rbind(data_frame(beta = betas[15],
thetas = thetas,
ll = llmat[15,])) %>%
rbind(data_frame(beta = betas[50],
thetas = thetas,
ll = llmat[50,])) %>%
rbind(data_frame(beta = betas[75],
thetas = thetas,
ll = llmat[75,])) %>%
rbind(data_frame(beta = betas[100],
thetas = thetas,
ll = llmat[100,])) %>%
mutate(beta = factor(beta)) %>%
ggplot(aes(x=thetas, y= ll, color = beta )) + geom_line() + theme_bw() +
geom_vline(aes(xintercept = .5),linetype = 2, color= 'black') +
ggtitle("Log - likelihood for possible values of theta, beta")data_frame(beta = betas,
theta = thetas[1],
ll = llmat[,1]) %>%
rbind(data_frame(beta = betas,
theta = thetas[25],
ll = llmat[,25])) %>%
rbind(data_frame(beta = betas,
theta = thetas[50],
ll = llmat[,50])) %>%
rbind(data_frame(beta = betas,
theta = thetas[100],
ll = llmat[,100])) %>%
mutate(theta = factor(theta)) %>%
ggplot(aes(x=beta, y= ll, color = theta )) + geom_line() + theme_bw() +
geom_vline(aes(xintercept = 1.2),linetype = 2, color= 'black') +
ggtitle("Log - likelihood for possible values of beta, theta")\[ p(\theta,\beta|Y,\alpha,\delta,\sigma) \propto p(y|...)p(\beta)p(\theta)\\ \beta \sim N(0,3)\\ \log(\theta) \sim N(1,1) \]
post <- map(1:length(betas),function(a) map_dbl(Xs,function(x) sum(dnorm(x = homog_subject_data$y,mean = Z %*% c(alpha,delta) + x$X_tilde * betas[a],
sd = sigma,log = F)) *
dnorm(x = betas[a],mean = 0,sd = 3,log = F) * 1)*
dlnorm(x = thetas[a],meanlog = 1, sdlog = 1, log =F))
postmat <- Reduce(f = rbind,x = post)par(mfrow=c(1,2))
persp3D(x = betas, y = thetas, z = postmat, phi = .6,theta = 0.5,xlab="beta",ylab="theta", zlab="density", ticktype = "detailed", main= "Posterior-unnormalized")
persp3D(x = betas, y = thetas, z = postmat, phi = 3, theta = 90, xlab="beta",ylab="theta", zlab="density", ticktype = "detailed")Looking at the \(\theta\) cross section, specifically
data_frame(beta = betas[1],
thetas = thetas,
post_dens = postmat[1,]) %>%
rbind(data_frame(beta = betas[50],
thetas = thetas,
post_dens = postmat[50,])) %>%
rbind(data_frame(beta = betas[100],
thetas = thetas,
post_dens = postmat[100,])) %>%
rbind(data_frame(beta = betas[120],
thetas = thetas,
post_dens = postmat[120,])) %>%
rbind(data_frame(beta = betas[150],
thetas = thetas,
post_dens = postmat[150,])) %>%
mutate(beta = factor(beta)) %>%
ggplot(aes(x=thetas, y= post_dens, color = beta )) + geom_line() + theme_bw() +
geom_vline(aes(xintercept = .5),linetype = 2, color= 'black') +
ggtitle("Posterior Density (unnormalized) for possible values of theta, beta")