RのMCMC修正案

Question

私は集団遺伝学のためにMCMCに取り組んできましたが、疑問があります。 私は統計に慣れていません。そのために私は困難があります。

MCMCを実行するコードが1000回あります。まず、0（50個の列= 50個の個体と1000個の繰り返しの1000行）の行列を作成します。 次に、ランダムなベクトルを作成して、行列の最初の行を置き換えます。このベクターには、1および2があり、集団1または集団2を表します。 また、遺伝子型頻度と50人の遺伝子型があります。 私が望むのは、遺伝子型の頻度と遺伝子型に応じて、個体がどの個体に属するのかを決定することです。 次に、ランダムな個体に割り当てられた母集団を変更し、新しい値を受け入れるかどうかを確認します。

niter <- 1000 
z <- matrix(0,nrow=niter,ncol=ncol(targetinds)) 
z[1,] <- sample(1:2, size=ncol(z), replace=T) 
lhood <- numeric(niter) 
lhood[1] <- compute_lhood_K2(targetinds, z[1,], freqPops) 
accepted <- 0 
priorz <- c(1e-6, 0.999999) 

for(i in 2:niter) { 

    z[i,] <- z[i-1,] 

    # propose new vector z, by selecting a random individual, proposing a new zi value 
    selind <- sample(1:nind, size=1) 
    # proposal probability of selecting individual at random 
    proposal_ratio_ind <- log(1/nind)-log(1/nind) 

    # propose a new index for the selected individual 
    if(z[i,selind]==1) { 
     z[i,selind] <- 2 
    } else { 
     z[i,selind] <- 1 
    } 

    # proposal probability of changing the index of individual is 1/2 
    proposal_ratio_cluster <- log(1/2)-log(1/2) 
    propratio <- proposal_ratio_ind+proposal_ratio_cluster 

    # compute f(x_i|z_i*, p) 
    # the probability of the selected individual given the two clusters 
    probindcluster <- compute_lhood_ind_K2(targetinds[,selind],freqPops) 
    # likelihood ratio f(x_i|z_i*,p)/f(x_i|z_i, p) 
    lhoodratio <- probindcluster[z[i,selind]]-probindcluster[z[i-1,selind]] 

    # prior ratio pi(z_i*)/pi(z_i) 
    priorratio <- log(priorz[z[i,selind]])-log(priorz[z[i-1,selind]]) 

    # accept new value according to the MH ratio 
    mh <- lhoodratio+propratio+priorratio 
    # reject if the random value is larger than the MH ratio 
    if(runif(1)>exp(mh)) { 
     z[i,] <- z[i-1,] # keep the same z 
     lhood[i] <- lhood[i-1] # keep the same likelihood 
    } else { # if accepted 
     lhood[i] <- lhood[i-1]+lhoodratio # update the likelihood 
     accepted <- accepted+1 # increase the number of accepted 
    } 
} 

新たに提案された値が尤度に比例するように、提案確率を変更する必要があるかどうか尋ねられます。これは、おそらくギブスのサンプリングMCMCアルゴリズムにつながると思われる。

私はこれを行うためにコードで何を変更するか分かりません。私はまた、提案の確率の概念と事前の選択方法を非常によく理解していません。

疑問を明確にする方法を知っている人に感謝します。

Answer 1

あなたの現在の提案はここで行われます。

# propose a new index for the selected individual 
    if(z[i,selind]==1) { 
     z[i,selind] <- 2 
    } else { 
     z[i,selind] <- 1 
    } 

個人がクラスタ1に割り当てられている場合は、クラスタ2（またはその逆）に割り当てることで確定的割り当てを切り替えることを提案します。

あなたは何であるかfreqPopsたちを示さなかったが、あなたはfreqPopsに応じて提案したいならば、私は上記のコードは、（少なくとも、それは私がときに理解するものである

z[i,selind] <- sample(c(1,2),size=1,prob=freqPops) 

で交換しなければならないと考えています尤度に基づいて提案したいと言いますが、あなたの声明は不明です）。 @papgeo、あなたの答えのための

proposal_ratio_cluster <- log(freqPops[z[i-1,selind]])-log(fregPops[z[i,selind]])