Approximate Survivor Function \(R_A(t)\)

The approximate survivor function is defined, and used, for \(t > 2\tau\):

\[\mathcal{R}_A(u) \backsim \sum_{i=1}^{k_A}\frac{\vec{c}_i\vec{r}_i}{\vec{r}_iW'(s_i)\vec{c}_i} e^{u s_i}\]

where \(c_i\) and \(r_i\) are the column and row eigenvectors of a \(H(s_i)\) (defined below), \(s_i\) are the roots of \(\mathop{det}W(s)\).

\[\begin{split}W(s) = sI - H(s)\\\end{split}\]
\[H(s) = \mathcal{Q}_{AA} + \mathcal{Q}_{AF}\left[sI-\mathcal{Q}_{FF}\right]^{-1} \left(I-e^{-(sI-\mathcal{Q}_{FF})\tau}\right)\mathcal{Q}_{FA}\]
\[W'(s) = I+\mathcal{Q}_{AF}\left[\bar{S}_{FF}(s)\left(sI-\mathcal{Q}_{FF}\right)^{-1} -\tau \left(I -\bar{\mathbf{S}}_{FF}(s)\right)\right]\bar{\mathcal{G}}_{FA}(s)\]
\[\bar{\mathbf{S}}_{FF}(s) = \int_0^\tau\!\mathrm{d}\,t\ e^{-st}e^{\mathcal{Q}_{FF}t}\]
\[\bar{\mathcal{G}}_{FA}(s) = \left[sI-\mathcal{Q}_{FF}\right]^{-1}\mathcal{Q}_{FA}\]

The approximate survivor function can be initialized from a \(\mathcal{Q}\)-matrix and the resolution \(\tau\):

python:
from dcprogs.likelihood import QMatrix, ApproxSurvivor

# 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

survivor = ApproxSurvivor(qmatrix, 1e-4);

print(survivor)
c++11:
#include <iostream>

#include <likelihood/approx_survivor.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::ApproxSurvivor survivor(qmatrix, 1e-4);

  std::cout << survivor << std::endl;

  return 0;
}

Note

The initialization shown above hides many details of the implementation, most notably how the roots are obtained. However, this object can also be initialed in exactly the same way as MissedEventsG, with the same sets parameters (except for nmax).

The open and shut time survivor likelihood can be computed using a single call:

python:

The python bindings accept both scalars and array inputs.

print("AF values\n"  \
      "---------\n")
times = [1e-4, 1e-5, 2.0e-5, 2.5e-5]
af_values = survivor.af(times)
for t, v in zip(times, af_values): print("  * at time t={0}:\n{1}\n".format(t, v))
c++11:
  std::cout << "AF values\n"
               "---------\n\n";
  std::cout << "  * at time t=" << 1e-4 <<":\n    "
            << DCProgs::numpy_io(survivor.af(1e-4), "    ") << "\n"
            << "  * at time t=" << 1.5e-4 <<":\n    " 
            << DCProgs::numpy_io(survivor.af(1.5e-4), "    ") << "\n"
            << "  * at time t=" << 2.0e-4 <<":\n    " 
            << DCProgs::numpy_io(survivor.af(2.0e-4), "    ") << "\n"
            << "  * at time t=" << 2.5e-4 <<":\n    "
            << DCProgs::numpy_io(survivor.af(2.5e-4), "    ") << "\n\n";

The coefficient and the exponents of the exponentials that make up the asymptotic expression are exposed as shown below.

python:
print("  * Exponents: {0}".format([root for matrix, root in survivor.af_components]))
c++11:
  std::cout << " * Exponents: ";
  for(DCProgs::t_uint i(0); i < survivor.nb_af_components(); ++i) 
    std::cout << std::get<1>(survivor.get_af_components(i)) << " ";
  std::cout << std::endl;