Exact Survivor Function \(R_A(t)\)

The exact survivor function is defined as:

\[\mathcal{R}_A(t) = \sum_{m=0}^\infty (-1)^m M_m(t-m\tau)\]

with,

\[\begin{split}M_m(t)=\left\{\begin{split} 0&\quad\text{if}\quad t < 0\\ \sum_{i=0}^k B_{im}(t) e^{-\lambda_i t}&\quad\text{if}\quad t \geq 0 & \end{split}\right.,\\ B_{im}(t) = \sum_{r=0}^mC_{imr}t^r,\end{split}\]

where the matrices \(C_{imr}\) are defined as a recursion:

\[\begin{split}\left\{\begin{eqnarray} C_{imm} &=& \frac{1}{m}D_iC_{i(m-1)(m-1)}\\ C_{im(l\leq m-1)} &=& \frac{1}{l}D_iC_{i(m-1)(l-1)} -\sum_{j\neq i}^k \sum_{r=l}^{m-1}\frac{r!}{l!(\lambda_i-\lambda_j)^{r-l+1}}D_jC_{i(m-1)r}\\ C_{im0} &=& \sum_{j\neq i}^k\sum_{r=0}^{m-1} \frac{r!}{(\lambda_i-\lambda_j)^{r+1}} \left[D_iC_{j(m-1)r}-D_jC_{i(m-1)r}\right] \end{eqnarray}\right.\end{split}\]

The initial values are \(C_{i00} = A_{iAA}\), with \(A_i\) the spectral decomposition of the \(\mathcal{Q}\)-matrix, and \(\lambda_i\) are its eigenvalues. Finally, the matrices \(D_i\) are defined as:

\[D_i = A_{iAF}e^{\mathcal{Q}_{FF}\tau}\mathcal{Q}_{FA}\]

Note

This recursion is implemented in the header likelihood/recursion.h in such a way that it can act upon a variety of objects. This makes testing it somewhat easier, since we can defined the \(D_i\), for instance, as scalars rather than matrices.

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

Python:

from HJCFIT.likelihood import QMatrix, ExactSurvivor

# 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 = ExactSurvivor(qmatrix, 1e-4);

print(survivor)

C++11:

#include <iostream>

#include <likelihood/exact_survivor.h>
#include <likelihood/errors.h>
 
int main() {

  // Define parameters.
  HJCFIT::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;
  HJCFIT::QMatrix qmatrix(matrix, /*nopen=*/2);

  HJCFIT::ExactSurvivor survivor(qmatrix, 1e-4);

  std::cout << survivor << std::endl;
                  << HJCFIT::numpy_io(survivor.recursion_af(i, m, l), "    ") << "\n\n";

  return 0;
}

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\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    "
            << HJCFIT::numpy_io(survivor.af(1e-4), "    ") << "\n"
            << "  * at time t=" << 1.5e-4 <<":\n    " 
            << HJCFIT::numpy_io(survivor.af(1.5e-4), "    ") << "\n"
            << "  * at time t=" << 2.0e-4 <<":\n    " 
            << HJCFIT::numpy_io(survivor.af(2.0e-4), "    ") << "\n"
            << "  * at time t=" << 2.5e-4 <<":\n    "

The details of the recursions, i.e. the \(C_{iml}\) matrices, can be accessed directly as shown below.

Python:

print("AF recusion matrices\n"  \
      "--------------------\n\n")
for i in range(5):
  for m in range(1, 3):
    for l in range(0, m+1):
      print( "  * C_{{{0}{1}{2}}}:\n{3}\n" \
             .format(i, m, l, survivor.recursion_af(i, m, l)) )

C++11:

  std::cout << "AF recusion matrices\n"
               "--------------------\n\n";
  for(HJCFIT::t_uint i(0); i < 5; ++i)
    for(HJCFIT::t_uint m(1); m < 3; ++m)
      for(HJCFIT::t_uint l(0); l <= m; ++l) 
        std::cout << "  * C_{" << i << m << l << "}:\n    "