MCMC en R Modifier la proposition

Question

J'ai travaillé avec MCMC pour la génétique des populations et j'ai quelques doutes. Je n'ai pas l'expérience des statistiques et à cause de cela j'ai de la difficulté.

J'ai le code pour exécuter MCMC, 1000 itérations. Je commence par créer une matrice avec 0 (50 colonnes = 50 individus et 1000 lignes pour 1000 itérations). Ensuite, je crée un vecteur aléatoire pour substituer la première ligne de la matrice. Ce vecteur a 1 et 2, représentant la population 1 ou la population 2. J'ai aussi des fréquences génotypiques et les génotypes des 50 individus. Ce que je veux c'est, selon les fréquences génotypiques et les génotypes, déterminer à quelle population appartient un individu. Ensuite, je vais continuer à changer la population assignée à un individu aléatoire et vérifier si la nouvelle valeur devrait être acceptée.

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

Il est demandé que je dois changer la probabilité de proposition afin que les nouvelles valeurs proposées sont proportionnelles à la probabilité. Cela conduit à un algorithme MCMC d'échantillonnage de Gibbs, soi-disant.

Je ne sais pas quoi changer dans le code pour le faire. Je ne comprends pas très bien le concept de probabilité de proposition et comment choisir le précédent. Reconnaissant si quelqu'un sait comment clarifier mes doutes.

Answer 1

Votre proposition actuelle se fait ici:

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

si l'individu est affecté au groupe 1, vous proposez de changer l'affectation déterministe en les affectant à la classe 2 (et vice versa).

Vous ne nous avez pas montrer ce que freqPops est, mais si vous voulez proposer selon freqPops alors je crois que le code ci-dessus doit être remplacé par

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

(au moins c'est ce que je comprends quand vous dites que vous voulez proposer en fonction de la probabilité - cependant, votre déclaration n'est pas claire).

Pour cela maintenant être un algorithme gibbs de MCCM valide d'échantillonnage que vous devez également modifier la ligne de code suivante:

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