シミュレーションの流れ

母集団の種類数$K$、$\mathbf{F}_{K}$を与えて、サンプルをランダムに作る。できたサンプルから母集団を推定し、$\hat{\mathbf{F}}_{K}$がどのくらいの場所にあるかを調べる。

基準にCDFを使うことにすると、$F_{K}$の分布はベータ分布になっているのでpbetaで簡単に調べられる。((気付いていなかったがその1の分布もdbeta(seq(0, 1,length.out=100) ,U[i], sum(U)-U[i]-k-2)でもっと簡単に計算できた))

上の議論により、$p$は $$ \begin{align} p &= &\int_{0}^{F} \text{pdf}_{\text{total}} \text{df}\\ &= &\displaystyle \sum_{K_{\text{sampled}} \leq K \lt \infty} \text{Cdf}_{K}(F) \end{align} $$

で定義され、pbetaを利用して計算する。

K.ori <- 4
Fs.ori <- c(0.5, 0.3, 0.1, 0.1)
n.sample <- 100
f.sample <- 100
K.sup <- 100
n.sim <- 2000

for (i in 1:n.sim){
  # サンプル作成
  U.sample <- as.vector(rmultinom(1, n.sample, Fs.ori))
  U.sample <- c(U.sample, 0)
  Ks <- seq(K.ori, K.sup)
  
  # ベータ分布の係数
  a <- U.sample + 1
  b <- outer(Ks, sum(U.sample) - U.sample - 1)
  
  # Kごとの分位点の計算
  p.Ks <- matrix(NA, length(Fs.ori), length(b[,1]))
  for(i in 1:length(Fs.ori)){
    p.Ks[i,] <- pbeta(Fs.ori[i], a[i], t(b)[i,])
    p.Ks[i,] <- ifelse(is.na(p.Ks[i,]), 0, p.Ks[i,])
  }
#    matplot(t(p.Ks), type='l', lty=1)

  # Kの尤度で重み付けした合計(ret)
  pmf.Ks <- exp(lgamma(sum(U.sample) + Ks) - lgamma(2 * sum(U.sample) + Ks))
  pmf.Ks <- pmf.Ks / sum(pmf.Ks)
#   matplot(t(t(p.Ks) * pmf.Ks), type='l', lty=1)
  ret <- apply(t(t(p.Ks) * pmf.Ks), 1, sum)

  # 書き出し
   cat(ret, '\n', file=paste('c(0.5, 0.3, 0.1, 0.1)', '-', 'Kmax=', K.sup, '.data', sep=''), append=T)
}

でサンプルを集められる。

評価

以前書いたプログラムで(51, 28, 11, 10)のサンプルが得られたときにpの値を計算すると、 0.508317, 0.699528, 0.340909, 0.464538であり、今回のプログラムでは0.54021207, 0.71517421, 0.34641109, 0.47016987が得られた。(==> その3 へ)

tak0kadaの何でもノート

発声練習、生存確認用。

種類数の予測(その2) - 「thinkbayes」第15章

シミュレーションの流れ

評価