rubias (version 0.3.4)

gsi_em_1: EM algorithm from the simplest GSI model for pi and the individual posterior probabilities

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)

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

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.

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)