gsi_em_1: EM algorithm from the simplest GSI model for pi and the...

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

Description

Using a matrix of scaled likelihoods, this function does an EM algorithm to climb the likelihood surface for pi, and computes the plug-in estimate of the posteriors for all the individuals. It returns the output in a list.

Usage

gsi_em_1(SL, Pi_init, max_iterations, tolerance, return_progression)

Arguments

SL	a 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 collections).
Pi_init	Starting value for the pi (collection mixture proportion) vector.
max_iterations	the maximum total number of reps iterations to do.
tolerance	the EM-algorithm will be considered converged when the sum over the elements of pi of the absolute value of the difference between the previous and the current estimate is less than tolerance.
return_progression	If true, then the pi_trace component of the output shows the value of pi visited en route to the end.

Value

gsi_em_1 returns a final Maximum-Likelihood estimate for pi and PofZ, as well as the number of iterations needed to reach convergence ("iterations_performed"), and traces of the pi values and change in pi in each iteration

Examples

# this is shown with a scaled likelihood matrix from self-assignment
# of the reference individuals

# 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))
test_em <- gsi_em_1(SL,
                    rep(1/params$C, params$C),
                    max_iterations = 10^6,
                    tolerance = 10^-7,
                    return_progression = TRUE)

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