Open mitzimorris opened 5 years ago
was file stan-dev/stan/src/docs/users-guide/R/funnel-fit.R:
funnel_fit <- stan(file='../../../models/misc/funnel/funnel_reparam.stan')
funnel_samples <- extract(funnel_fit,permuted=TRUE,inc_warmup=FALSE);
funnel_df <- data.frame(x1=funnel_samples$x[,1],y=funnel_samples$y[])
ggp <-
ggplot(funnel_df,aes(x=x1,y=y)) +
coord_cartesian(xlim=c(-20,20), ylim=c(-9,9)) +
scale_x_continuous("x[1]", expand=c(0,0), breaks=c(-20,-10,0,10,20)) +
scale_y_continuous("y", expand=c(0,0), breaks=c(-9,-6,-3,0,3,6,9)) +
labs(title="Funnel Samples (transformed model)\n") +
geom_point(shape='.', alpha=1, color="black")
png(filename="funnel-fit.png", width=1200,height=1200,res=300);
print(ggp);
dev.off();
was file stan-dev/stan/src/docs/users-guide/R/funnel-plot.R:
library("ggplot2");
p_funnel <- function(x1,y) {
return(dnorm(y,0.0,3.0,log=TRUE) + dnorm(x1,0.0,exp(y/2),log=TRUE));
}
K = 200;
x1 <- rep(0,(K+1)^2);
y <- rep(0,(K+1)^2);
log_p_y_x1 <- rep(0,(K+1)^2);
pos <- 1;
for (m in 1:(K+1)) {
for (n in 1:(K+1)) {
y[pos] <- -9.0 + 18.0 * (m - 1) / K;
x1[pos] <- -20.0 + 40.0 * (n - 1) / K;
log_p_y_x1[pos] <- p_funnel(x1[pos],y[pos]);
pos <- pos + 1;
}
}
library("scales");
df <- data.frame(x1=x1,
y=y,
log_p_y_x1=log_p_y_x1);
funnel_plot <-
ggplot(df, aes(x1,y,fill = log_p_y_x1)) +
labs(title = "Funnel Density (log scale)\n") +
geom_tile() +
scale_x_continuous("x[1]",expand=c(0,0),limits=c(-20,20), breaks=c(-20,-10,0,10,20)) +
scale_y_continuous(expand=c(0,0),limits=c(-9,9), breaks=c(-9,-6,-3,0,3,6,9)) +
scale_fill_gradient2("log p(y,x[1])\n",
limits=c(-18,2),
midpoint=-6.75, mid="lightyellow",
low="gray95", high="darkblue", na.value="transparent",
breaks=c(0,-8,-16),
labels=c("0","-8","-16"));
png(filename="funnel.png", width=1500,height=1200,res=300);
print(funnel_plot);
dev.off();
was stan-dev/stan/src/docs/users-guide/R/non-identified-plot.R:
library(ggplot2);
library(grid);
N <- 100;
y <- rnorm(N);
loc_sum <- function(lambda1,lambda2) {
sum(dnorm(y,lambda1 + lambda2, 1, log=TRUE));
}
loc_sum_prior <- function(lambda1,lambda2) {
sum(dnorm(y,lambda1 + lambda2, 1, log=TRUE)) +
sum(dnorm(c(lambda1,lambda2),0,1,log=TRUE));
}
one_param <- function(xs) {
dta <- rnorm(100,0,1);
result <- rep(NA, length(xs));
i <- 1;
for (x in xs) {
result[i] <- sum(dnorm(dta,xs[i],1,log=TRUE));
i <- i + 1;
}
return(result);
}
p_one_param <-
ggplot(data.frame(mu=c(-20, 20)), aes(mu)) +
labs(title = "Proper Posterior (without Prior)\n") +
stat_function(fun=one_param) +
labs(x=expression(mu), y="log p") +
theme(aspect.ratio=1,
panel.border=element_blank(),
plot.margin=unit(c(0,0,0,0),"cm"),
text=element_text(size=28),
axis.title=element_text(size=32));
png(res=100,height=800,width=900,file="one_param_identified.png");
print(p_one_param);
dev.off();
K <- 500;
ub <- 25;
lambda_1 <- seq(-ub,ub,len=K);
lambda_2 <- seq(-ub,ub,len=K);
v_lambda_1 <- rep(NA,K^2);
v_lambda_2 <- rep(NA,K^2);
v_density <- rep(NA,K^2);
# use prior (loc_sum_prior)
pos <- 1;
for (m in 1:K) {
for (n in 1:K) {
v_lambda_1[pos] <- lambda_1[m];
v_lambda_2[pos] <- lambda_2[n];
v_density[pos] <- loc_sum_prior(lambda_1[m], lambda_2[n]);
pos <- pos + 1;
}
}
df_vals <- list(lambda_1=v_lambda_1,lambda_2=v_lambda_2,log_density=v_density);
df_prior <- as.data.frame(df_vals);
# plot and save as png
p_prior <-
ggplot(df_prior,aes(x=lambda_1,y=lambda_2,fill=log_density)) +
labs(title = "Proper Posterior (with Prior)\n") +
geom_tile() +
labs(x=expression(lambda[1]), y=expression(lambda[2])) +
scale_fill_gradient2("log p",
limits=c(-600,-120),
midpoint=-400, low="gray95", high="darkblue",
mid="lightyellow", na.value="transparent",
breaks=c(-200, -400, -600),
labels=c("-200", "-400", "-600")) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
theme(aspect.ratio=1,
panel.border=element_blank(),
plot.margin=unit(c(0,0,0,0),"cm"),
text=element_text(size=28),
axis.title=element_text(size=32),
legend.title=element_text(size=30,face="plain",color="gray25"),
legend.text=element_text(size=20,color="gray40"),
legend.margin=unit(0.3,"cm") )
png(res=100,height=800,width=900,file="non-identified-plus-prior.png");
plot(p_prior)
dev.off();
plot(p_prior)
# no prior (loc_sum)
pos <- 1;
for (m in 1:K) {
for (n in 1:K) {
v_density[pos] <- loc_sum(lambda_1[m], lambda_2[n]);
pos <- pos + 1;
}
}
df_vals <- list(lambda_1=v_lambda_1,lambda_2=v_lambda_2,log_density=v_density);
df_no_prior <- as.data.frame(df_vals);
# plot and save as png
p_no_prior <-
ggplot(df_no_prior,aes(x=lambda_1,y=lambda_2,fill=log_density)) +
labs(title = "Improper Posterior (without Prior)\n") +
geom_tile() +
labs(x=expression(lambda[1]), y=expression(lambda[2])) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
scale_fill_gradient2("log p",
limits=c(-600,-120),
midpoint=-400, low="gray95", high="darkblue",
mid="lightyellow", na.value="transparent",
breaks=c(-200, -400, -600),
labels=c("-200", "-400", "-600") ) +
theme(aspect.ratio=1,
panel.border=element_blank(),
plot.margin=unit(c(0,0,0,0),"cm"),
text=element_text(size=28),
axis.title=element_text(size=32),
legend.title=element_text(size=30,face="plain",color="gray25"),
legend.text=element_text(size=20,color="gray40"),
legend.margin=unit(0.3,"cm") )
theme(aspect.ratio=1,text=element_text(size=28),
axis.title=element_text(size=32),
legend.margin=unit(0.3,"cm"));
png(res=100,height=800,width=900,file="non-identified.png")
print(p_no_prior)
dev.off();
Summary:
markdown code includes png (speeds up build), but we don't have the R code checked in which generates these images.
Description:
following files include pngs - need code which creates these
Additional Information:
Provide any additional information here.
Current Version:
v2.18.0