Ideal Likelihood \(\mathcal{G}(t)\)¶
A wrapper around QMatrix is provided which allows the calculation of the ideal likelihood:
\[\begin{split}\mathcal{G}(t) = \left(
\begin{eqnarray}
e^{-t\mathcal{Q}_{af}} & 0 \\
0& e^{-t\mathcal{Q}_{fa}} \\
\end{eqnarray}
\right)\end{split}\]
This object can be initialized directly from a QMatrix.
| python: | from dcprogs.likelihood import QMatrix, IdealG, expm
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)
idealG = IdealG(qmatrix)
print(idealG)
|
|---|---|
| c++11: | #include <iostream>
#include <exception>
#include <likelihood/idealG.h>
#include <likelihood/errors.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::IdealG idealG(qmatrix);
std::cout << idealG << std::endl;
return 0;
}
|
It provides the ideal likelihood as a function of time, as well as the laplace transforms:
| python: | from numpy import dot, identity, abs, all
idealG_fa = dot(expm(2e-4 * qmatrix.ff), qmatrix.fa)
assert all( abs(idealG_fa - idealG.fa(2e-4)) < 1e-8 )
inversion = -0.5 * identity(2) - qmatrix.aa
assert all( abs( dot(inversion, idealG.laplace_af(-0.5)) - qmatrix.af ) < 1e-8 )
Note This package can be compiled to use real number of 128bits. However, despite providing arrays for reals of that size, numpy does not include the corresponding linear algebra functions. A few useful functions, such as expm in this example, are provided to remediate to this situation. |
|---|---|
| c++11: | DCProgs::t_rmatrix const idealG_fa = (2e-4*qmatrix.ff()).exp()*qmatrix.fa();
if( ((idealG.fa(2e-4) - idealG_fa).array().abs() > 1e-8).any() )
throw DCProgs::errors::Runtime("Not so ideal idealG");
DCProgs::t_rmatrix const inversion = -0.5 * DCProgs::t_rmatrix::Identity(2, 2) - qmatrix.aa();
if( ((inversion * idealG.laplace_af(-0.5) - qmatrix.af()).array().abs() > 1e-8).any() )
throw DCProgs::errors::Runtime("Not so ideal idealG");
|