Maximum Subinterval Likelihood
The maximum likelihood estimator MIL (Qin et al 1996) is extended to handle data recorded with a varying stimulus. We consider the stimulus as piece-wise constant, and adjust the transition matrices each time it changes. This may split some events, so we derive a transition probability matrix that allows subsequent events in the same class. Thus the new name, maximum subinterval likelihood.
See the source code: browse
MIL Background
In MIL, the conditional likelihood of a dwell in class a and the following transition to class b is given as
$$\begin{equation} {}_e\!G_{ab}(t) = e^{{}_e\!Q_{aa}t'}{}_e\!Q_{ab} \label{eq:eGab} \end{equation}$$where
$$\begin{equation} t' = t - t_d \label{eq:tp} \end{equation}$$ $$\begin{equation} {}_e\!Q_{aa} = Q_{aa} - Q_{a\bar{a}}(I - e^{Q_{\bar{a}\bar{a}}t_d})Q_{\bar{a}\bar{a}}^{-1}Q_{\bar{a}a} \label{eq:eQaa} \end{equation}$$ $$\begin{equation} {}_e\!Q_{ab} = e^{t_dQ_{a\bar{a}}(I - e^{Q_{\bar{a}\bar{a}}t_d})Q_{\bar{a}\bar{a}}^{-1}Q_{\bar{a}a}} \times [Q_{ab} - Q_{ac}(I - e^{Q_{cc}t_d})Q_{cc}^{-1}Q_{cb}] \label{eq:eQab} \end{equation}$$\({}_e\!Q_{ab}\) is calculated differently in modern QUB:
$$\begin{equation} {}_e\!Q_{ab} = [Q_{ab} - Q_{ac}(I - e^{Q_{cc}t_d})Q_{cc}^{-1}Q_{cb}]e^{Q_{bb}t_d} \label{eq:eQab2} \end{equation}$$The forward function is
$$\begin{equation} \alpha^{\tau}_k = \alpha^{\tau}_{k-1}{}_e\!G_{a_ka_{k+1}}(t_k) \label{eq:alpha} \end{equation}$$and
$$\begin{equation} LL = log(\alpha^{\tau}_L 1) \label{eq:ll} \end{equation}$$Subinterval Likelihood
Following the lead of Colquhoun and Hawkes (1982), eq \eqref{eq:eGab} combines the conditional probability of staying in class a (at least apparently) with the probability per second of going to class b. The resulting probabilities per second are often greater than 1, so their use as likelihoods in MIL leads to surreal "LL". This is mainly just an annoyance, but it does prevent use of Aikaike and Bayes information criteria.
A more serious difficulty for MSL is that the stimulus may change mid-dwell, and with it the Q matrix. We need a formula that allows for subsequent dwells of the same class. In equation \eqref{eq:eGab}, \(e^{{}_e\!Q_{aa}t'}\) does not depend on the event ending after time t. So the likelihood of a partial dwell is
$$\begin{equation} L_a(t) = e^{{}_e\!Q_{aa}(t - t_d)} \label{eq:La} \end{equation}$$The final \(t_d\) of each event is spent mysteriously in transition, whether back to the same class or into a new one. This is the conditional likelihood of being in class b after exactly \(t_d\):
$$\begin{equation} Adead = e^{Q t_d} \label{eq:Adead} \end{equation}$$ $$\begin{equation} L_{ab} = Adead_{ab} = Adead_{i \in a,j \in b} \label{eq:Lab} \end{equation}$$Put them together in a new forward function and get LL from \eqref{eq:alpha}:
$$\begin{equation} \alpha^{\tau}_k = \alpha^{\tau}_{k-1} L_a(t_k) L_{a_ka_{k+1}} \label{eq:alphatau} \end{equation}$$Constrained Parameters
This is the algorithm used by MIL, Mac etc to transform between k, a column vector of \(N_k\) rate constants, constrained by
$$\begin{equation} A_{in} \cdot k = B_{in} \label{eq:Ain} \end{equation}$$and the free parameters \(p \in \Re^{N_p}\) used by the optimizer. A and B are derived so that
$$\begin{equation} k = A \cdot p + B, and \label{eq:k} \end{equation}$$ $$\begin{equation} p = A^T \cdot k \label{eq:p} \end{equation}$$In the final section I extend the method to scale all free parameters to start at 1.
Clean up Constraints
Remove linearly dependent constraints
The Singular Value Decomposition (SVD) technique below doesn't work well when constraints duplicate each other; that is, none of the augmented vectors \(c_i^T = \left [ A_{ij} \vdots B_i \right ]\) may be a linear combination of the others. To find and eliminate these constraints, build C out of rows \(c_i\): $$\begin{equation} C = \left [ A \vdots B \right ]^T \label{eq:c} \end{equation}$$
define the partitions:
$$\begin{equation} C_{i..j} = \left [ {}^{c_{1i}}\ddots_{c_{N_k+1 j}} \right ]; C_{i..i} = c_i \label{eq:Cij} \end{equation}$$and with i from 2 increasing, remove any column \(c_j\), where \(i \lt j\), if
$$\begin{equation} C_{1..i} \cdot c'_j = c_j \label{eq:C1i} \end{equation}$$can be solved for \(c'_j\). This can be tested with SVD:
$$\begin{equation} U \cdot W \cdot V^T = SVD(C_{1..i}) \label{eq:svd} \end{equation}$$ $$\begin{equation} C^-1 \approx V \cdot W^{-1} \cdot U^T \label{eq:Cinv} \end{equation}$$* if \(W_{ii} = 0\), define \(W_{ii}^{-1} = 0\)
If \(c_j = c''_j\), then \(c_j\) is the linear combination \(c'_j\) of the previous constraints, so remove it.
Alternately, or as a double-check, the \(c_i\) are linearly independent if \(C_{1..N_c}\) has rank \(N_c\). A constraint \(c_j\), where \(j = i + 1\) can be removed if \(rank(C_{1..i}) = rank(C_{1..j})\).
Again with SVD: the singular values of \(C_{1..i}\) are the nonzero diagonal entries of W in \eqref{eq:svd}. \(rank(C_{1..i})\) is the number of singular values.
Look for Inconsistent Constraints
Consider \(rank(A_{in}^T) = rank(A_{in})\). If \(rank(A_{in}) \lt rank(C)\), there is a constraint that can be expressed as a linear combination of others, except with a different right-hand-side (\(B_{in}\)) component, a contradiction. The user must resolve this. We could identify which constraints conflict with each other using a procedure similar to the last section.
Solve the Linear System
Now we have \(A_{in}\), \(B_{in}\) with \(N_c\) independent and compatible constraints. To solve the system \eqref{eq:Ain} for the subspace of constrained k, observe that any such k can be written
$$\begin{equation} k = k_h + k_p \label{eq:khp} \end{equation}$$where \(k_h\) is in the null space of \(A_{in}\); i.e. \(A_{in} \cdot k_h = 0\), and \(k_p\) is a particular vector satisfying \eqref{eq:Ain}.
The null space
Using SVD,
$$\begin{equation} U \cdot W \cdot V^T = A_{in} \label{eq:svdnull} \end{equation}$$Now, assuming SVD sorts the singular values non-increasing, the columns of
$$\begin{equation} A = \left [ V_{R+1 \cdots N_k} \right ] \label{eq:Avr} \end{equation}$$form a basis for the null space of \(A_{in}\), so
$$\begin{equation} k_h = A \cdot p = 0 \text{ for all } p \in \Re^{N_k-R} \label{eq:kh} \end{equation}$$The Pseudo-Inverse
Find the particular nonhomogeneous solution \(k_p\) by inverting \eqref{eq:Ain} using SVD pseudo-inversion:
$$\begin{equation} A_{in} \cdot k_p = B_{in} \label{eq:Ainkp} \end{equation}$$ $$\begin{equation} A_{in}^{-1} \approx V W^{-1} U^T \label{eq:Aininv} \end{equation}$$ $$\begin{equation} A_{in}^{-1} A_{in} k_p \approx V W^{-1} U^T B_{in} \label{eq:Aininvkp} \end{equation}$$ $$\begin{equation} B = k_p \approx V W^{-1} U^T B_{in} \label{eq:B} \end{equation}$$Putting it together,
$$\begin{equation} k = A \cdot p + B \label{eq:kdup} \end{equation}$$ $$\begin{equation} \frac{\partial k}{\partial p} = A \label{eq:dkdp} \end{equation}$$Scale All Parameters to 1
The DFP optimizer works poorly when components of p have different magnitude. To scale them so they're all ones at the starting rate constants \(k0\), define the vector \(s\) of scaled parameters
$$\begin{equation} p = A_s \cdot s + B_s \label{eq:pAs} \end{equation}$$ $$\begin{equation} s = A^{-1}_s (p - B_s) \label{eq:sAinv} \end{equation}$$ $$\begin{equation} \frac{\partial p}{\partial s} = S \label{eq:dpds} \end{equation}$$where \(A_s\) is the \(N_p x N_p\) diagonal matrix with entries
$$\begin{equation} p0 = A^T \cdot k0 \label{eq:p0At} \end{equation}$$ $$\begin{equation} A_{s,ii} = \begin{cases} 1.0, & \mbox{if }p0_i\mbox{ is 0} \\ p0_i, & \mbox{otherwise} \end{cases} \label{eq:Asii} \end{equation}$$and \(B_s\) is the column vector with entries