複数（> 2）色のポリアールプロセスをモデル化するにはどうすればよいですか？

Question

私はPolya Urnモデルを構築しようとしています。 2色はうまく行きましたが、3色では問題が発生しました。

ndraws<-1000; nexps<-2000; Distribution.yellow<-matrix(0,ndraws,1); for (k in  1:nexps){ 
red<- 1; 
yellow<- 1; 
blue<-1 ; 
for (n in 1:ndraws){ 
    drawn<-sample(0:2,size=1,prob=c(red,yellow,blue)/(red+yellow +blue)) 
    red<-?? ; 
    blue<-?? ; 
    yellow<-?? ; 
    } 
    Distribution.yellow[k]<-yellow/(red+yellow+blue) } 

私の問題は、このコード行を変換することにある：骨壷に追加し、それぞれの余分なボールの中へ

drawn<-sample(0:2,size=1,prob=c(red,yellow,blue)/(red+yellow +blue)) 

を。 （したがって、疑問符）。

drawn<-sample(0:1,size=1,prob=c(red,blue)/(red+blue)) 
red<-red+(1-drawn); 
blue<-blue+(drawn); 

しかし、以上の2色があるとき、これは明らかに動作しません：次のように私はそれをやった2色で

。どのように私は3つ以上の色でアプローチする必要がありますか？ Wikipediaによれば

Answer 1

、ポリA URN処理するためのルールは、次のとおり

1つのボールが骨壷からランダムに引き出され、その色が観察され、それを壷に戻し、同じ色の追加の玉を壷に加え、選択プロセスを繰り返す。

つまり、ボールを描くと、そのコロ（u）rのボールの数が1ずつ増えます。

だから我々は1つの赤いボールであればdrawn==0、1黄色のボールであればdrawn==1、それ以外の場合は一つの青のボールを追加if文を...セットアップすることができ

ndraws <- 1000 
nexps <- 500 
set.seed(101) 
yellow_final <- numeric(nexps) 
for (k in 1:nexps) { 
    red <- 1; yellow <- 1; blue<-1 
    for (n in 1:ndraws) { 
     drawn <- sample(0:2,size=1,prob=c(red,yellow,blue)/(red+yellow+blue)) 
     if (drawn==0) { 
      red <- red+1 
     } else if (drawn==1) { 
      yellow <- yellow+1 
     } else blue <- blue+1 
    } 
    yellow_final[k]<-yellow/(red+yellow+blue) 
} 

写真：

par(las=1,bty="l") 
hist(yellow_final,col="gray",freq=FALSE, 
    xlab="Prop. yellow after 1000 draws") 

Answer 2

一般化された溶液：

ndraws<-1000 
nexps<-2000 
colors <- c('red', 'blue', 'yellow') # add balls with other colors 
initial.num.balls <- c(1,1,1) # can have different numbers of balls to start with 
ball.to.observe <- 'yellow' 
distribution.ball.to.observe <- replicate(nexps, { 
    urn <- rep(colors, initial.num.balls) # polya's urn 
    count.balls <- as.list(initial.num.balls) 
    names(count.balls) <- colors 
    for (i in 1:ndraws) { 
    drawn <- sample(urn, 1) 
    count.balls[[drawn]] <- count.balls[[drawn]] + 1 
    urn <- c(urn, drawn) 
    } 
    count.balls[[ball.to.observe]]/sum(as.numeric(count.balls)) 
}) 
library(ggplot2) 
ggplot() + stat_density(aes(distribution.ball.to.observe), bw=0.01)