3

ベイジアン更新の理解を助けるために、bayesianbiologistのコードを使用してきました。アニメーション プロットの作成方法を学ばなければならないので、更新のアニメーション プロットを作成するのは楽しい練習になると思いました。これは予想以上に厳しいものでした。この問題に関する 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() を挿入しようとしましたが、これはプロットの追加的な性質と競合すると思います。

これらの方法の一方または両方を適切に機能させる方法を知っている人はいますか? これは私にとってはちょっとしたおもちゃの例ですが、プロットをアニメーション化する方法を理解することが役立つ/必要な場合に、プロットする必要がある興味深いアニメーションがいくつかあります。

助けてくれてありがとう！

4

2 に答える 2

1

MatthewK が言ったように、あなたは への呼び出しを置き忘れましたpdf。ただし、これで起動して実行できるようになります。

sim_bayes <- function(p=0.5, N=10, y_lim=15, prior_a=1, prior_b=1) {
    success <- 0
    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

        # Start a new plot.
        curve(dbeta(x,prior_a,prior_b), lty=2,
              xlim=c(0,1), ylim=c(0,y_lim), xlab='p', ylab='Posterior Density')
        # Update plot.
        curve(dbeta(x,success+prior_a, (i-success) + prior_b), add=TRUE)

        legend('topright',
               legend=c('Prior','Updated Posteriors','Final Posterior'),
               lty=c(2,1,1), col=c('black','black','red'))

        dev.off()
    }
}

# `x` had no visible binding in your implementation, so I took the following
# from the `dbeta` documentation example.
x <- seq(0, 1, length=21)
sim_bayes()
于 2012-08-19T08:37:23.737 に答える
1

ここで作成するプロットは 1 つだけです。add=TRUE現在のプロットに追加する場合は、新しいプロットを作成しないでください。したがって、それを削除するとうまくいくはずです：

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)) #plot updated
      }
curve(dbeta(x,success+prior_a,(i-success)+prior_b),col='red',lwd=1.5) #plot final posterior
}

テスト： 

sim_bayes(p=0.6,N=90,prior_a=1,prior_b=1)

複数のプロットを提供します

于 2012-08-19T08:44:05.280 に答える