CODEWITHSUNDEEP
rmcmc
Targaryel
Targaryel

Reputation: 117

MCMC in R Modify Proposal

I've been working with MCMC for population genetics and I have some doubts. I'm not experienced in statistics and because of that I have difficulty.

I have code to run MCMC, 1000 iterations. I start by creating a matrix with 0's (50 columns = 50 individuals and 1000 lines for 1000 iterations). Then I create a random vector to substitute the first line of the matrix. This vector has 1's and 2's, representing population 1 or population 2. I also have genotype frequencies and the genotypes of the 50 individuals. What I want is to, according to the genotype frequencies and genotypes, determine to what population an individual belongs. Then, I'll keep changing the population assigned to a random individual and checking if the new value should be accepted.

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

It is asked that I have to change the proposal probability so that the new proposed values are proportional to the likelihood. This leads to a Gibbs sampling MCMC algorithm, supposedly.

I don't know what to change in the code to do this. I also don't understand very well the concept of proposal probability and how to chose the prior.

Grateful if someone knows how to clarify my doubts.

Upvotes: 0

Views: 120

Answers (1)

papgeo
papgeo

Reputation: 473

Your current proposal is done here:

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

if the individual is assigned to cluster 1, then you propose to switch assignment deterministically by assigning them to cluster 2 (and vice versa).

You didn't show us what freqPops is, but if you want to propose according to freqPops then I believe the above code has to be replaced by 

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

(at least that is what I understand when you say you want to propose based on the likelihood - however, that statement of yours is unclear).

For this now to be a valid mcmc gibbs sampling algorithm you also need to change the next line of code:

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

Upvotes: 2

Related Questions