# (C) Petteri Pajunen 2006 # # # bgauss.R # # demonstrates Bayesian inference on mean value # # of a Normal distribution with known variance. X11(width=5,height=5) t <- seq(-10,10,by=0.05) sigma <- 1 # known standard deviation mu <- 5 # unknown mean value xreal <- dnorm(t,mean=mu,sd=sigma) # data distribution # h=plot(t,xreal); plot(t,xreal,ylim=c(0,1),lwd=2,type="l",col="red",main="TRUE DATA DISTRIBUTION") # xlabel('RED=TRUE DISTRIBUTION'); cat("Choose the prior mean: ") mu0 <- scan("",n=1) cat("Choose the prior std: ") sigma0 <- scan("",n=1) xprior <- dnorm(t,mean=mu0,sd=sigma0) lines(t,xprior,lwd=2,col="blue") # xlabel('RED=TRUE PDF, BLUE=PRIOR PDF'); # initialize posterior pdf with prior pdf sigmap <- sigma0 mup <- mu0 for (i in 1:5) { y <- rnorm(1,mean=mu,sd=sigma) points(y,0.2,col="red",lwd=3) # h=plot(y,0.2,'rs','MarkerSize',10,'MarkerFaceColor','r'); # calculate new posterior pdf using data y newprecision <- sigmap^(-2)+sigma^(-2) mup <- (sigmap^(-2)*mup+sigma^(-2)*y)/(newprecision) sigmap <- newprecision^(-1/2) xpost <- dnorm(t,mean=mup,sd=sigmap) plot(t,xreal,type="l",ylim=c(0,1),lwd=2,col="red",main="RED=TRUE, BLUE=PRIOR, GREEN=POSTERIOR") lines(t,xprior,lwd=2,col="blue") points(y,0.2,col="red",lwd=3) lines(t,xpost,lwd=2,col="green") # xlabel('RED=TRUE PDF, BLUE=PRIOR PDF, GREEN=POSTERIOR PDF'); #cat('Press any key to see next sample') scan(quiet=TRUE) #if (i<5) # delete(h2); #end; } X11(width=5,height=5) mupred <- mup sigmapred <- sqrt(sigma^2+sigmap^2) xpred <- dnorm(t,mean=mupred,sd=sigmapred) plot(t,xreal,type="l",ylim=c(0,1),lwd=2,col="red",main="RED=TRUE, BLUE=PRIOR, GREEN=PREDICTIVE") lines(t,xprior,lwd=2,col="blue") points(y,0.2,col="red",lwd=3) # lines(t,xpost,lwd=2,col="green") lines(t,xpred,lwd=2,col="green") # xlabel('RED=TRUE PDF, BLUE=PRIOR PDF, GREEN=POSTERIOR PDF, DASHED=PRED PDF'); # fprintf('CORRECT MEAN: 5, CORRECT STD: 1\n'); # fprintf('POSTERIOR MEAN: %.2f, POSTERIOR STD: %.2f\n',mup,sigmap); scan(quiet=TRUE) dev.off() dev.off()