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:
- from the knowledge of the roots and the determinant equations
- directly from a \(Q\)-matrix, using the default root-finding mechanism
- from a \(Q\)-matrix, using a custom root-finding mechanism (c++ only)
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 \({}^eG_{AF}(t)\) and \({}^eG_{FA}(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\tau\)). 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.