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

In rubias: Bayesian Inference from the Conditional Genetic Stock Identification Model

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)

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.

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.

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)

rubias documentation built on Feb. 10, 2022, 1:06 a.m.