Likelihood \(L(Q)\)

The callable object \(L_{\log10}(Q)\) provides an interface designed to compute the log10-likelihood as a function of the rates, from the knowledge of the observed open and shut time intervals and the (type of) initial and final occupancies. This object is especially suited for maximizing the likelihood of a given mechanism. It is defined as

\[L_{\log10}(Q) = \frac{1}{\mathrm{ln} 10}\sum_b \mathrm{ln} L_b(Q)\]

with \(L_b(Q)\) defined here, the index \(b\) refers to the individual bursts, and \(\mathrm{ln}\) is the natural logarithm.

The purpose of this class is to provide an interface for maximizing the likelihood. It computes, for a given set of observed open and shut intervals, the likelihood \(\frac{\ln(L(Q))}{ln 10}\), where \(L(Q)\) is declared in the likelihood equation.

A callable object \(L(Q)\) exists in both c++ and python. It can be initialized as follows

python:
from dcprogs.likelihood import Log10Likelihood

bursts = [  [0.1, 0.2, 0.1],                  # 1st burst 
            [0.2],                            # 2nd burst
            [0.15, 0.16, 0.18, 0.05, 0.1] ]   # 3rd burst
""" List of bursts.  

    Each burst is a list of observed open and shut intervals. 
    There should always be an odd number of intervals, since bursts end in a shut states.
"""

likelihood = Log10Likelihood(bursts, nopen=2, tau=0.01, tcritical=None)

print(likelihood)
c++11:
#include <iostream>
#include <likelihood/likelihood.h>

int main() {

  DCProgs::t_Bursts bursts{
     {0.1, 0.2, 0.1},                  /* 1st burst */
     {0.2},                            /* 2nd burst */
     {0.15, 0.16, 0.18, 0.05, 0.1}     /* 3rd burst */
  };
  
  DCProgs::Log10Likelihood likelihood( bursts, /*nopen=*/2, /*tau=*/1e-2,
                                       /*tcrit=*/DCProgs::quiet_nan );
 
  std::cout << likelihood << std::endl;
  return 0;
}

The initialization of bursts above is done in using two newer c++11 coding techniques: initializer lists, and uniform initialization. It may not be available from all compilers just yet...

Once the objects are initialized, the input attributes can be accessed (and modified) through the ‘.’ operator: likelihood.nopen = 2.

Note

Log10Likelihood() uses equilibrium occupancies depending on the value of its attribute tcritical:

  • if it is None, numpy.NaN, or negative, then the equilibrium occupancies are used
  • if it a strictly positive real number, then the CHS vectors are computed

Similarly, in c++, tcritical can be set to DCProgs::quiet_nan to trigger calculations with equilibrium occupancies.

It is required that the bursts have been pre-filtered so that there are no intervals smaller than the resolution \(\tau\). This can be done using time_filter() (python) or interval_filter() (python).

The likelihood for any Q-matrix can then be computed by calling the likelihood objects as though they were function. The following snippets are inserted at the tail end of the previous code.

python:
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 ] ]

result = likelihood(matrix)

print("Computation: {0}".format(result))

The function can take any sort square matrix, whether using standard python lists or a numpy array. It can only take one matrix at a time.
c++11:
  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::t_real const result = likelihood(matrix);

The return is the log-likelihood associated with the bursts and the input Q-matrix. In both python and c++, the functions accepts either a matrix or an actual DCProgs::QMatrix (python) object. In the former case, the number of open states is set to nopen.

It should be noted that the python the bursts are accessed in python directly from the likelihood using normal sequence operations. Only a small subset of sequence operations where implemented.

python:
from numpy import all, abs, NaN
# Get the number of bursts
assert len(likelihood) == 3
# Check the second burst is what we expect
assert all(abs(likelihood[2] - [0.15, 0.16, 0.18, 0.05, 0.1]) < 1e-12)
# Modify, and check the second burst
likelihood[2] = [0.15, 0.16, 0.18]
assert all(abs(likelihood[2] - [0.15, 0.16, 0.18]) < 1e-12)
# Add an extra burst, and check that it is there
likelihood.append([0.25, 0.013, 0.013])
assert len(likelihood) == 4
assert all(abs(likelihood[3] - [0.25, 0.013, 0.013]) < 1e-12)
c++11:

DCProgs::Log10Likelihood::bursts is a public member and can be accessed directly.

Finally, some of the attributes, namely, Log10Likelihood.tcritical, Log10Likelihood.upper_bound, Log10Likelihood.lower_bound, act both as parameters and as switch when given special values. These special values are None and numpy.NaN in python and DCProgs::quiet_nan in c++. In python, the special values will always be transformed to None.

python:
assert likelihood.tcritical is None
likelihood.tcritical = NaN
assert likelihood.tcritical is None
likelihood.tcritical = -1
assert likelihood.tcritical is None
likelihood.tcritical = 0.5
assert likelihood.tcritical is not None and abs(likelihood.tcritical - 0.5) < 1e-12

assert likelihood.lower_bound is None
likelihood.lower_bound = NaN
assert likelihood.lower_bound is None
likelihood.lower_bound = -1e6
assert likelihood.lower_bound is not None and abs(likelihood.lower_bound + 1e6) < 1e-12
likelihood.lower_bound = None
assert likelihood.lower_bound is None