Processing math: 100%

DCProgs 0.9 documentation

Exact Survivor Function RA(t)

«  Occupancies   ::   Contents   ::   Approximate Survivor Function RA(t)  »

Exact Survivor Function RA(t)

The exact survivor function is defined as:

RA(t)=m=0(1)mMm(tmτ)

with,

Mm(t)={0ift<0ki=0Bim(t)eλitift0,Bim(t)=mr=0Cimrtr,

where the matrices Cimr are defined as a recursion:

{Cimm=1mDiCi(m1)(m1)Cim(lm1)=1lDiCi(m1)(l1)kjim1r=lr!l!(λiλj)rl+1DjCi(m1)rCim0=kjim1r=0r!(λiλj)r+1[DiCj(m1)rDjCi(m1)r]

The initial values are Ci00=AiAA, with Ai the spectral decomposition of the Q-matrix, and λi are its eigenvalues. Finally, the matrices Di are defined as:

Di=AiAFeQFFτQFA

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 Di, for instance, as scalars rather than matrices.

The survivor function can be initialized from a Q-matrix and the resolution τ:

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;
}

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

«  Occupancies   ::   Contents   ::   Approximate Survivor Function RA(t)  »