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:

python:	from dcprogs.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. 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::ExactSurvivor survivor(qmatrix, 1e-4); std::cout << survivor << std::endl; << DCProgs::numpy_io(survivor.recursion_af(i, m, l), " ") << "\n\n"; return 0; }

from dcprogs.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.
  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::ExactSurvivor survivor(qmatrix, 1e-4);

  std::cout << survivor << std::endl;
                  << DCProgs::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:

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 " << 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 "

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    "
            << 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    "

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

python:

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(DCProgs::t_uint i(0); i < 5; ++i) for(DCProgs::t_uint m(1); m < 3; ++m) for(DCProgs::t_uint l(0); l <= m; ++l) std::cout << " * C_{" << i << m << l << "}:\n "

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(DCProgs::t_uint i(0); i < 5; ++i)
    for(DCProgs::t_uint m(1); m < 3; ++m)
      for(DCProgs::t_uint l(0); l <= m; ++l) 
        std::cout << "  * C_{" << i << m << l << "}:\n    "

DCProgs 0.9 documentation

Exact Survivor Function \(R_A(t)\)

Exact Survivor Function \(R_A(t)\)¶