treetime/clock_tree.py
def get_max_posterior_region(self, node, fraction = 0.9):
        '''
        If temporal reconstruction was done using the marginal ML mode, the entire distribution of
        times is available. This function determines the interval around the highest
        posterior probability region that contains the specified fraction of the probability mass.
        In absense of marginal reconstruction, it will return uncertainty based on rate
        variation. If both are present, the wider interval will be returned.

        Parameters
        ----------

         node : PhyloTree.Clade
            The node for which the posterior region is to be calculated

         interval : float
            Float specifying who much of the posterior probability is
            to be contained in the region

        Returns
        -------
         max_posterior_region : numpy array
            Array with two numerical dates delineating the high posterior region

        '''
        if node.marginal_inverse_cdf=="delta":
            return np.array([node.numdate, node.numdate])


        min_max = (node.marginal_pos_LH.xmin, node.marginal_pos_LH.xmax)
        min_date, max_date = [self.date2dist.to_numdate(x) for x in min_max][::-1]
        if node.marginal_pos_LH.peak_pos == min_max[0]: #peak on the left
            return self.get_confidence_interval(node, (0, fraction))
        elif node.marginal_pos_LH.peak_pos == min_max[1]: #peak on the right
            return self.get_confidence_interval(node, (1.0-fraction ,1.0))
        else: # peak in the center of the distribution
            rate_contribution = self.date_uncertainty_due_to_rate(node, ((1-fraction)*0.5, 1.0-(1.0-fraction)*0.5))

            # construct height to position interpolators left and right of the peak
            # this assumes there is only one peak --- might fail in odd cases
            from scipy.interpolate import interp1d
            from scipy.optimize import minimize_scalar as minimize
            pidx = np.argmin(node.marginal_pos_LH.y)
            pval = np.min(node.marginal_pos_LH.y)
            left =  interp1d(node.marginal_pos_LH.y[:(pidx+1)]-pval, node.marginal_pos_LH.x[:(pidx+1)],
                            kind='linear', fill_value=min_max[0], bounds_error=False)
            right = interp1d(node.marginal_pos_LH.y[pidx:]-pval, node.marginal_pos_LH.x[pidx:],
                            kind='linear', fill_value=min_max[1], bounds_error=False)

            # function to minimize -- squared difference between prob mass and desired fracion
            def func(x, thres):
                interval = np.array([left(x), right(x)]).squeeze()
                return (thres - np.diff(node.marginal_cdf(np.array(interval))))**2

            # minimza and determine success
            sol = minimize(func, bracket=[0,10], args=(fraction,))
            if sol['success']:
                mutation_contribution = self.date2dist.to_numdate(np.array([right(sol['x']), left(sol['x'])]).squeeze())
            else: # on failure, return standard confidence interval
                mutation_contribution = None

            return self.combine_confidence(node.numdate, (min_date, max_date),
                                      c1=rate_contribution, c2=mutation_contribution)
def get_max_posterior_region(self, node, fraction = 0.9):
        '''
        If temporal reconstruction was done using the marginal ML mode, the entire distribution of
        times is available. This function determines the interval around the highest
        posterior probability region that contains the specified fraction of the probability mass.
        In absense of marginal reconstruction, it will return uncertainty based on rate
        variation. If both are present, the wider interval will be returned.

        Parameters
        ----------

         node : PhyloTree.Clade
            The node for which the posterior region is to be calculated

         interval : float
            Float specifying who much of the posterior probability is
            to be contained in the region

        Returns
        -------
         max_posterior_region : numpy array
            Array with two numerical dates delineating the high posterior region

        '''
        if node.marginal_inverse_cdf=="delta":
            return np.array([node.numdate, node.numdate])


        min_max = (node.marginal_pos_LH.xmin, node.marginal_pos_LH.xmax)
        min_date, max_date = [self.date2dist.to_numdate(x) for x in min_max][::-1]
        if node.marginal_pos_LH.peak_pos == min_max[0]: #peak on the left
            return self.get_confidence_interval(node, (0, fraction))
        elif node.marginal_pos_LH.peak_pos == min_max[1]: #peak on the right
            return self.get_confidence_interval(node, (1.0-fraction ,1.0))
        else: # peak in the center of the distribution
            rate_contribution = self.date_uncertainty_due_to_rate(node, ((1-fraction)*0.5, 1.0-(1.0-fraction)*0.5))

            # construct height to position interpolators left and right of the peak
            # this assumes there is only one peak --- might fail in odd cases
            from scipy.interpolate import interp1d
            from scipy.optimize import minimize_scalar as minimize
            pidx = np.argmin(node.marginal_pos_LH.y)
            pval = np.min(node.marginal_pos_LH.y)
            left =  interp1d(node.marginal_pos_LH.y[:(pidx+1)]-pval, node.marginal_pos_LH.x[:(pidx+1)],
                            kind='linear', fill_value=min_max[0], bounds_error=False)
            right = interp1d(node.marginal_pos_LH.y[pidx:]-pval, node.marginal_pos_LH.x[pidx:],
                            kind='linear', fill_value=min_max[1], bounds_error=False)

            # function to minimize -- squared difference between prob mass and desired fracion
            def func(x, thres):
                interval = np.array([left(x), right(x)]).squeeze()
                return (thres - np.diff(node.marginal_cdf(np.array(interval))))**2

            # minimza and determine success
            sol = minimize(func, bracket=[0,10], args=(fraction,))
            if sol['success']:
                mutation_contribution = self.date2dist.to_numdate(np.array([right(sol['x']), left(sol['x'])]).squeeze())
            else: # on failure, return standard confidence interval
                mutation_contribution = None

            return self.combine_confidence(node.numdate, (min_date, max_date),
                                      c1=rate_contribution, c2=mutation_contribution)
If temporal reconstruction was done using the marginal ML mode, the entire distribution of
        times is available. This function determines the interval around the highest
        posterior probability region that contains the specified fraction of the probability mass.
        In absense of marginal reconstruction, it will return uncertainty based on rate
        variation. If both are present, the wider interval will be returned.

        Parameters
        ----------

         node : PhyloTree.Clade
            The node for which the posterior region is to be calculated

         interval : float
            Float specifying who much of the posterior probability is
            to be contained in the region

        Returns
        -------
         max_posterior_region : numpy array
            Array with two numerical dates delineating the high posterior region