Given a lineage tree and distribution information, calculate the likelihood of the tree and choose analysis of highest likelihood

[shak360]

The following code uses LL calculations for parameter estimation. We can easily write a LL for the entire tree assuming no transitions.

    def bernoulliParameterEstimatorNumerical(self):
        '''Estimates the Bernoulli parameter for a given population using MLE numerically'''
        population = self.group # assign population to a variable
        fate_holder = [] # instantiates list to hold cell fates as 1s or 0s
        for cell in population: # go through every cell in the population
            if not cell.isUnfinished(): # if the cell has lived a meaningful life and matters
                fate_holder.append(cell.fate*1) # append 1 for dividing, and 0 for dying

        def negLogLikelihoodBern(locBernGuess, fate_holder):
            """ Calculates the log likelihood for bernoulli. """
            return -1*np.sum(sp.bernoulli.logpmf(k=fate_holder, p=locBernGuess))

        res = minimize(negLogLikelihoodBern, x0=0.5, bounds=((0,1),), method="SLSQP", args=(fate_holder))

        return res.x[0]

    def gompertzParameterEstimatorNumerical(self):
        '''Estimates the Gompertz parameters for a given population using MLE numerically'''
        population = self.group # assign population to a variable
        tau_holder = [] # instantiates list
        for cell in population: # go through every cell in the population
            if not cell.isUnfinished(): # if the cell has lived a meaningful life and matters
                tau_holder.append(cell.tau) # append the cell lifetime

        def negLogLikelihoodGomp(gompParams, tau_holder):
            """ Calculates the log likelihood for gompertz. """
            return -1*np.sum(sp.gompertz.logpdf(x=tau_holder,c=gompParams[0], scale=gompParams[1]))

        res = minimize(negLogLikelihoodGomp, x0=[2,50], bounds=((0,5),(0,None)), method="SLSQP", options={'maxiter': 1e7}, args=(tau_holder))

return res.x

[shak360]

Sorry this likelihood calculation is without transition rates.

[shak360]

See #197