# Manually computing a posterior ################################## rm(list=ls()) set.seed(2391) posterior<-function(theta, N=50, k=20){ like<-(theta^(k))*(1-theta)^(N-k) pr<-ifelse(abs(theta)<0.2, 1, 0) out<-pr*like out<-ifelse(theta>1, 0, out) out<-ifelse(theta<0, 0, out) return(out) } theta.seq<-seq(from=0, to=0.4, by=0.001) post.den<-posterior(theta.seq) post.den.plot<-(post.den/sum(post.den))*length(theta.seq) plot(post.den.plot~theta.seq, type="l") # Cromwell's Rule: no prior density = no posterior density # recompute, but with a "flat" (and improper) prior density rm(list=ls()) set.seed(2391) posterior<-function(theta, N=50, k=20){ like<-(theta^(k))*(1-theta)^(N-k) pr<-5 out<-pr*like out<-ifelse(theta>1, 0, out) out<-ifelse(theta<0, 0, out) return(out) } theta.seq<-seq(from=0, to=0.8, by=0.001) post.den<-posterior(theta.seq) post.den.plot<-post.den*(length(theta.seq)/sum(post.den)) plot(post.den.plot~theta.seq, type="l") # recompute with an "informative" prior density rm(list=ls()) set.seed(2391) posterior<-function(theta, N=50, k=20){ like<-(theta^(k))*(1-theta)^(N-k) pr<-dnorm(theta, mean=0.3, sd=0.05) out<-pr*like out<-ifelse(theta>1, 0, out) out<-ifelse(theta<0, 0, out) return(out) } theta.seq<-seq(from=0, to=0.8, by=0.001) post.den<-posterior(theta.seq) post.den.plot<-post.den*(length(theta.seq)/sum(post.den)) plot(post.den.plot~theta.seq, type="l") lines(dnorm(theta.seq, mean=0.3, sd=0.05)~theta.seq, lty=2) # investigate the beta prior density theta.seq<-seq(from=0, to=1, by=0.001) plot(dbeta(theta.seq, shape1=10, shape2=30)~theta.seq, type="l") lines(dbeta(theta.seq, shape1=3, shape2=3)~theta.seq, lty=2) # recompute with a beta prior density rm(list=ls()) set.seed(2391) posterior<-function(theta, N=50, k=20, alpha=10, beta=30){ like<-(theta^(k))*(1-theta)^(N-k) pr<-(theta^alpha)*(1-theta)^(beta) out<-pr*like out<-ifelse(theta>1, 0, out) out<-ifelse(theta<0, 0, out) return(out) } theta.seq<-seq(from=0, to=0.8, by=0.001) post.den<-posterior(theta.seq) post.den.plot<-post.den*(length(theta.seq)/sum(post.den)) plot(post.den.plot~theta.seq, type="l") lines(dbeta(theta.seq, shape1=10, shape2=30)~theta.seq, lty=2) # bayes factors ################################# 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) z<-runif(100) model.1<-glm(y~1) model.2<-glm(y~x) model.3<-glm(y~x+z) exp(as.numeric(logLik(model.2)-logLik(model.1))) exp(as.numeric(logLik(model.3)-logLik(model.2)))