MCMC в предложении изменения R

Я работаю с MCMC в области популяционной генетики, и у меня есть некоторые сомнения. У меня нет опыта в статистике, и из-за этого у меня возникают трудности.

У меня есть код для запуска MCMC, 1000 итераций. Я начинаю с создания матрицы с нулями (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 выборки Гиббса.

Я не знаю, что для этого изменить в коде. Я также не очень хорошо понимаю концепцию вероятности предложения и то, как выбрать априорную.

Благодарен, если кто-то знает, как развеять мои сомнения.


r mcmc
person Targaryel    schedule 23.06.2017    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Ваше текущее предложение сделано здесь:

# 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]])
person papgeo    schedule 24.06.2017
comment
Большое спасибо за ответ, @papgeo! Что я действительно хочу сделать, так это предложить в соответствии со значениями правдоподобия, рассчитанными для отдельного человека. 2 значения правдоподобия для этого человека, по одному для каждой популяции. С вашим ответом я думаю, что смогу сделать эту часть, я просто использовал значения правдоподобия в prob, чтобы наиболее вероятное было выбрано больше раз. - person Targaryel; 24.06.2017
comment
Что касается второй части, не должно ли fregPops [z [i, selind]] быть первым, а freqPops [z [i-1, selind]] - вторым? Итак, это: log (freqPops [z [i, selind]]) - log (fregPops [z [i-1, selind]])? Кроме того, не могли бы вы объяснить мне или дать объяснение Гиббса в подобном случае? Простите, что так много спрашиваю! - person Targaryel; 24.06.2017
comment
Я считаю, что порядок, который я предложил, правильный. Речь идет о исправлении возможного дисбаланса в предложении. Так что изменение предложения не должно сильно повлиять на результат. Но изменение приоритета может иметь эффект. На данный момент настоятель говорит, что первая группа маловероятна. 1e-6. - person papgeo; 24.06.2017