Как смоделировать несколько (> 2) цветных урнов?

Question

Я пытаюсь построить модель Поля Урна. Два цвета пошли хорошо, с тремя, однако я столкнулся с некоторыми неприятностями.

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); 

Но это, очевидно, не работает, когда есть более двух цветов. Как мне подойти с тремя или более цветами?

Answer 1

Согласно Wikipedia, правила для процесса урн Полна являются:

один мяч обращается случайным образом из урны и ее цвет наблюдается; затем он возвращается в урну, а в урну добавляется дополнительный шар того же цвета, и процесс выбора повторяется.

Другими словами, нанесение шарика увеличивает количество шаров этого колора (u) r на единицу.

Таким образом, мы можем создать if заявление, что добавляет один красный шар, если drawn==0, один желтый шар, если drawn==1 и один синий шар в противном случае ...

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)