# binomial learning ##################### rm(list=ls()) set.seed(324807) alpha = 2 # number of prior "successes" beta = 8 # number of prior "failures" x <- rbinom(n=50, size=1, prob=0.2) # the data post.samp <- rbeta(10000, shape1=(sum(x==1)+alpha-1), shape2=(sum(x==0)+beta-1)) # sample out of posterior hist(post.samp) theta.plot <- seq(from=0.1, to=0.4, by=0.001) # values of theta to plot post.plot <- dbeta(theta.plot, shape1=(sum(x==1)+alpha-1), shape2=(sum(x==0)+beta-1) ) # posterior density values lines(200*post.plot ~ theta.plot, type="l", col="red") # compare to flat prior post.plot.2 <- dbeta(theta.plot, shape1=(sum(x==1)), shape2=(sum(x==0)) ) # posterior density values lines(200*post.plot.2 ~ theta.plot, type="l", col="blue") # mean estimate mean(post.samp) # variance/SEs can be calculated analytically... or: sqrt(var(post.samp)) # count learning ##################### rm(list=ls()) set.seed(6776) a <- 100 # prior sum of counts b <- 10 # prior sample size x <- rpois(n=200, lambda=6) # the data S <- sum(x) post.samp <- rgamma(10000, shape=S+a, scale=1/(b+200)) hist(post.samp) x<-seq(from=5.6, to=7, by=0.001) post.plot <- dgamma(x, shape=S+a, scale=1/(b+200)) lines(1000*post.plot ~ x, type="s", col="red") # compare to flat prior post.plot.2 <- dgamma(x, shape=S, scale=1/(200)) lines(1000*post.plot.2 ~ x, type="s", col="blue") mean(post.samp) sd(post.samp) # continuous/normal learning # known sigma ############################## rm(list=ls()) set.seed(6776) mu.0 <- 12 sigma.0 <- 0.5 n <- 200 x <- rnorm(n, mean=5, sd=1) mu.star <- ( (mu.0*sigma.0^(-2)) + mean(x)*(n/var(x)) ) / ( sigma.0^(-2) + n/(var(x)) ) sigma.star <-sqrt( (sigma.0^(-2) + n/(var(x)))^(-1) ) post.samp <- rnorm(10000, mean=mu.star, sd=sigma.star) hist(post.samp) x<-seq(from=4.8, to=5.5, by=0.001) post.plot <- dnorm(x, mean=mu.star, sd=sigma.star) lines(500*post.plot ~ x, type="s", col="red") x<-seq(from=4.8, to=5.5, by=0.001) post.plot <- dnorm(x, mean=mean(x), sd=sqrt(var(x)/(n-1))) lines(500*post.plot ~ x, type="s", col="blue") mean(post.samp) sd(post.samp) # bayesian regression analysis ################################# rm(list=ls()) set.seed(123921) # make a fake data set x<-runif(100) y<-0.5+2*x+rnorm(100, mean=0, sd=1) # bind the independent variables and the constant into a matrix X<-cbind(1,x) # run a model linear.model<-glm(y~x) # recover the information from this model coefs<-coefficients(linear.model) VCV<-vcov(linear.model) resid<-residuals(linear.model) # input prior information for beta coefficients coefs.prior<-c(0,1) b0<-coefs.prior VCV.prior<-matrix(data=c(0.75,0,0,0.75), ncol=2, nrow=2) B0<-VCV.prior # show the prior density for the beta coefficients library(mvtnorm) b.plot<-seq(from=-3, to=3, by=0.1) pdf.beta<-matrix(data=NA, nrow=length(b.plot), ncol=length(b.plot)) for(i in 1:length(b.plot)){ for(j in 1:length(b.plot)){ pdf.beta[i,j]<-dmvnorm(c(b.plot[i],b.plot[j]), mean=coefs.prior, sigma=VCV.prior) } } contour(x=b.plot, y=b.plot, z=pdf.beta) filled.contour(x=b.plot, y=b.plot, z=pdf.beta) filled.contour(x=b.plot, y=b.plot, z=pdf.beta, col=gray(21:0/21)) # input the prior information about the error process for the DGP (sigma) v0<-2 sigma0<-1 # show this prior on sigma library(MCMCpack) sigma.prior<-dinvgamma(seq(from=0, to=3, by=0.05), shape=v0/2, scale=v0*sigma0/2) plot(sigma.prior~(seq(from=0, to=3, by=0.05)), type="l") XtX<-t(X)%*%X b1<-qr.solve(qr.solve(B0)+XtX)%*%(qr.solve(B0)%*%b0+XtX%*%coefs) B1<-qr.solve(qr.solve(B0)+XtX) v1=v0+length(X) r<-t(b0-coefs)%*%qr.solve(B0+qr.solve(XtX))%*%(b0-coefs) v1s1<-v0*sigma0+sum(resid^2)+r # show the marginal posterior density for beta coefficient on x library(mnormt) b.plot<-seq(from=-6, to=6, by=0.1) posterior.beta<-matrix(data=NA, nrow=length(b.plot), ncol=length(b.plot)) for(i in 1:length(b.plot)){ for(j in 1:length(b.plot)){ posterior.beta[i,j]<-dmt(c(b.plot[i],b.plot[j]),mean=as.vector(b1), S=as.numeric(v1s1/v1)*B1, df=v1, log=T) } } filled.contour(x=b.plot, y=b.plot, z=posterior.beta, col=gray(30:0/30)) b.plot<-seq(from=0.5, to=3, by=0.01) x.beta<-dmt(b.plot,mean=as.vector(b1[2]), S=as.numeric(v1s1/v1)*B1[2,2], df=v1, log=F) plot(x.beta~b.plot, type="l") x.prior<-dnorm(b.plot,mean=as.vector(b0[2]), sd=sqrt(B0[2,2])) lines(x.prior~b.plot, lty=2)