rubias (version 0.3.4)

gsi_mcmc_1: MCMC from the simplest GSI model for pi and the individual posterior probabilities

Description

Using a matrix of scaled likelihoods, this function samples values of pi and the posteriors for all the individuals. It returns the output in a list.

Usage

gsi_mcmc_1(SL, Pi_init, lambda, reps, burn_in, sample_int_Pi, sample_int_PofZ)

Value

gsi_mcmc_1 returns a list of three. $mean lists the posterior means for collection proportions pi and for the individual posterior probabilities of assignment PofZ. $sd returns the posterior standard deviations for the same values.

If the corresponding sample_int variables are not 0, $trace contains samples taken from the Markov chain at intervals of sample_int_(variable) steps.

Arguments

SL

matrix of the scaled likelihoods. This is should have values for each individual in a column (going down in the rows are values for different populations).

Pi_init

Starting value for the pi (collection mixture proportion) vector.

lambda

the prior to be added to the collection allocations, in order to generate pseudo-count Dirichlet parameters for the simulation of a new pi vector

reps

total number of reps (sweeps) to do.

burn_in

how many reps to discard in the beginning when doing the mean calculation. They will still be returned in the traces if desired

sample_int_Pi

the number of reps between samples being taken for Pi traces. If 0 no trace samples are taken

sample_int_PofZ

the number of reps between samples being taken for the traces of posterior of each individual's origin. If 0 no trace samples are taken.

Examples

# this demonstrates it with scaled likelihoods computed from
# assignment of the reference samples

# we have to get the ploidies to pass to tcf2param_list
locnames <- names(alewife)[-(1:16)][c(TRUE, FALSE)]
ploidies <- rep(2, length(locnames))
names(ploidies) <- locnames

params <- tcf2param_list(alewife, 17, ploidies = ploidies)
logl <- geno_logL(params)
SL <- apply(exp(logl), 2, function(x) x/sum(x))
lambda <- rep(1/params$C, params$C)
mcmc <- gsi_mcmc_1(SL, lambda, lambda, 200, 50, 5, 5)