为了更好地理解贝叶斯更新,我一直在使用来自Bayesian生物学家的代码。由于我需要学会创建动画图形,所以我认为创建一个更新的动画图形会是一个有趣的练习。然而,这比我预期的要困难得多。受Rob Hyndman关于此问题的博客的启发,我尝试创建以下内容:
library(animation)
setwd("~/Dropbox/PriorUpdating") #set working directory
## Simulate Bayesian Binomial updating
sim_bayes<-function(p=0.5,N=10,y_lim=15,prior_a=1,prior_b=1)
{
print(paste("The prior expectation of p is ",prior_a/(prior_a+prior_b)))
success<-0
curve(dbeta(x,prior_a,prior_b),xlim=c(0,1),ylim=c(0,y_lim),xlab='p',ylab='Posterior Density',lty=2)
legend('topright',legend=c('Prior','Updated Posteriors','Final Posterior'),lty=c(2,1,1),col=c('black','black','red'))
for(i in 1:N)
{
if(runif(1,0,1)<=p) success<-success+1 #this is where we see if there is a "success"
curve(dbeta(x,success+prior_a,(i-success)+prior_b),add=TRUE) #plot updated
}
curve(dbeta(x,success+prior_a,(i-success)+prior_b),add=TRUE,col='red',lwd=1.5) #plot final posterior
}
oopt = ani.options(interval = 0)
saveMovie(sim_bayes(p=0.6,N=90,prior_a=1,prior_b=1),interval=0.1,width=580,height=400)
ani.options(oopt)
然而,这只生成了最终的图表。因此,我想尝试输出所有图表的PDF文件。
library(animation)
setwd("~/Dropbox/PriorUpdating") #set working directory
## Simulate Bayesian Binomial updating
sim_bayes<-function(p=0.5,N=10,y_lim=15,prior_a=1,prior_b=1)
{
print(paste("The prior expectation of p is ",prior_a/(prior_a+prior_b)))
success<-0
curve(dbeta(x,prior_a,prior_b),xlim=c(0,1),ylim=c(0,y_lim),xlab='p',ylab='Posterior Density',lty=2)
legend('topright',legend=c('Prior','Updated Posteriors','Final Posterior'),lty=c(2,1,1),col=c('black','black','red'))
for(i in 1:N)
{
pdf(paste("posterior",i,".pdf",sep=""),height=4,width=6.5)
if(runif(1,0,1)<=p) success<-success+1 #this is where we see if there is a "success"
curve(dbeta(x,success+prior_a,(i-success)+prior_b),add=TRUE) #plot updated
#print(paste(success,"successes and ",i-success," failures"))
dev.off()
}
pdf(paste("posterior_final",".pdf",sep=""),height=4,width=6.5)
curve(dbeta(x,success+prior_a,(i-success)+prior_b),add=TRUE,col='red',lwd=1.5) #plot final posterior
dev.off()
}
然而,这给我带来了以下错误。
Error in plot.xy(xy.coords(x, y), type = type, ...) :
plot.new has not been called yet
我尝试在某些位置插入plot.new(),但我认为这会与图的加性特性冲突。
有人有什么想法,如何使这两种方法都正常工作?虽然这对我来说是一个玩具示例,但我需要绘制一些更有趣的动画,而理解如何动画化图表将是有用/必要的。
感谢任何帮助!