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 算法。

我不知道要在代码中更改什么才能执行此操作。我也不是很明白proposal probability的概念以及如何选择先验

如果有人知道如何澄清我的疑惑，将不胜感激。

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 吉布斯采样算法，您还需要更改下一行代码：

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

R 中的 MCMC 修改提案

MCMC in R Modify Proposal

r

mcmc