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
    }
}

有人要求我必须更改提议概率，以便新提议的值与可能性成正比。据推测，这导致了 Gibbs 采样 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)

（至少当你说你想根据可能性提出建议时，我是这么理解的——但是，你的陈述不清楚）。

现在要成为有效的 mcmc gibbs 采样算法，您还需要更改下一行代码：

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