Processing math: 100%

DCProgs 0.9 documentation

Missed-Events Likelihood eG(t)

«  Likelihood L(Q)   ::   Contents   ::   Ideal Likelihood G(t)  »

Missed-Events Likelihood eG(t)

The callable object DCProgs::MissedEventsG provides an interface to compute the likelihood eG(t) of open and shut events as a function of their lengths, for a fixed Q-matrix. It has the ability to compute both exact and approximate missed-events likelihood, returning one or the other depending on a given time cutoff.

The asymptotic expression of the likelihood can be computed from the knowledge of the roots of a specific equations. On the one hand, root-finding can be a fairly difficult numerical operation. On the other, it would be more convenient if we can initialize DCProgs::MissedEventsG directly from a Q-matrix object. As such, there are several means to initialize the functor:

Initialization from a Q-matrix

python:
from numpy import all, abs, arange
from dcprogs.likelihood import QMatrix, DeterminantEq, MissedEventsG

# Define parameters.
qmatrix = QMatrix([ [-3050,        50,  3000,      0,    0], 
                    [2./3., -1502./3.,     0,    500,    0], 
                    [   15,         0, -2065,     50, 2000], 
                    [    0,     15000,  4000, -19000,    0], 
                    [    0,         0,    10,      0,  -10] ], 2)
tau = 1e-4
 
eG_automatic = MissedEventsG(qmatrix, tau)
c++11:
#include <iostream>
#include <exception>

#include <likelihood/missed_eventsG.h>
#include <likelihood/root_finder.h>
 
int main() {

  // Define parameters.
  DCProgs::t_rmatrix matrix(5 ,5);
  matrix << -3050,        50,  3000,      0,    0, 
            2./3., -1502./3.,     0,    500,    0,  
               15,         0, -2065,     50, 2000,  
                0,     15000,  4000, -19000,    0,  
                0,         0,    10,      0,  -10;
  DCProgs::QMatrix qmatrix(matrix, /*nopen=*/2);
  DCProgs::t_real const tau(1e-4); // in seconds

  DCProgs::MissedEventsG eG_automatic(qmatrix, tau);

  return 0;
}

A fair amount of work goes on behind the scene. First reasonable upper and lower bounds for the roots obtained (DCProgs::find_lower_bound_for_roots(), and DCProgs::find_upper_bound_for_roots()). Then intervals for each roots are computed (DCProgs::find_root_intervals()). And finally, the roots themselves are obtained (DCProgs::brentq()). All this work is done automatically in the case of this particular instantiation. A few extra parameters to control the root-finding process can be passed to the c++ and python constructors.

Initialization from the roots and determinant equations

python:
determinant_eq = DeterminantEq(qmatrix, tau)
af_roots = [( -3045.285776037674, 1), (-162.92946543451328, 1)]
fa_roots = [(-17090.192769236815, 1), (-2058.0812921673496, 1), (-0.24356535498785126, 1)]
eG_from_roots = MissedEventsG(determinant_eq, af_roots, determinant_eq.transpose(), fa_roots)
c++11:
  DCProgs::DeterminantEq determinant_eq(qmatrix, tau);
  std::vector<DCProgs::Root> af_roots{
    { /*root=*/ -3045.285776037674,   /*multiplicity=*/ 1}, 
    { /*root=*/ -162.92946543451328,  /*multiplicity=*/ 1}
  };
  std::vector<DCProgs::Root> fa_roots{
    { /*root=*/ -17090.192769236815,      /*multiplicity=*/ 1},
    { /*root=*/  -2058.0812921673496,     /*multiplicity=*/ 1},
    { /*root=*/     -0.24356535498785126, /*multiplicity=*/ 1}
  };
  DCProgs::MissedEventsG eG_from_roots( determinant_eq, af_roots, 
                                        determinant_eq.transpose(), fa_roots );

In this case, it is expected the roots are known, somehow, as well as their multiplicity. This format allows external root-finding methods to be interfaced with the packages.

Initialization from the Q-matrix and a root finding function

Given a root-finding function, it is possible to instantiate eG. The root finding function should take a determinant equation as input, and return a vector of DCProgs::Root as output. In the code below, we show how the prior initialization could be recreated.

  auto find_roots = [](DCProgs::DeterminantEq const &_det) {
    return DCProgs::find_roots(_det, 1e-12, 1e-12, 100, DCProgs::quiet_nan, DCProgs::quiet_nan);
  };
  DCProgs::MissedEventsG eG_from_func(qmatrix, tau, find_roots);

This is mostly a convenience function, to make it slightly easier to interface with other root-finding methods in c++. This interface is not explicitly available in python, although it can be created with ease.

Computing eGAF(t) and eGFA(t)

The computation of the likelihood matrices is illustrated by comparing the three initialization methods (two in python). It is clear that all three yield the same function.

python:
assert all(abs(eG_from_roots.af(tau) - eG_automatic.af(tau)) < 1e-8)
assert all(abs(eG_from_roots.fa(tau) - eG_automatic.fa(tau)) < 1e-8)

# Checks the three initialization are equivalent at different times
# The functions can be applied to arrays. 
times = arange(tau, 10*tau, 0.1*tau)
assert eG_from_roots.af(times).shape == (len(times), 2, 3)
assert eG_from_roots.fa(times).shape == (len(times), 3, 2)
assert all(abs(eG_from_roots.af(times) - eG_automatic.af(times)) < 1e-8)
assert all(abs(eG_from_roots.fa(times) - eG_automatic.fa(times)) < 1e-8)
c++11:
  for(DCProgs::t_real t(tau);  t < 10*tau; t += tau * 0.1) {

    if(    ((eG_from_roots.af(t) - eG_from_func.af(t)).array().abs() > 1e-8).any()  
        or ((eG_from_roots.fa(t) - eG_from_func.fa(t)).array().abs() > 1e-8).any() )
      throw DCProgs::errors::Runtime("root != func");

    if(    ((eG_from_roots.af(t) - eG_automatic.af(t)).array().abs() > 1e-8).any() 
        or ((eG_from_roots.fa(t) - eG_automatic.fa(t)).array().abs() > 1e-8).any() )
      throw DCProgs::errors::Runtime("root != automatic");
  }

The python bindings accept any input that can be transformed to a numpy array of reals. If the input is a scalar, then the AF and FA blocks are returned. If the input is an array, then an array of similar shape is returned, where each component is a matrix.

The DCProgs::MissedEventsG provides further functionality. For instance, the cutoff point between exact and asymptotic calculations can be set explicitly (it defaults to t<3τ). And the likelihood can be computed in Laplace space (see DCProgs::MissedEventsG::laplace_af and DCProgs::MissedEventsG::laplace_fa). We invite users to turn to the python and the c++ API for more details.

«  Likelihood L(Q)   ::   Contents   ::   Ideal Likelihood G(t)  »