View source for LChrisman/Computing Derivatives for MetaLog

This page outlines an algorithm for computing the <math>i^{th}</math> derivative of a Keelin MetaLog.

From Eq(6) in the Baucells, Chrisman, Keelin paper draft:
:<math>M^{(i)}(y) = \mu^{(i)}(y) + s^{(i)} \ell(y) + \sum_{j=0}^{i-1} s^{(j)}(y) \ell^{(i-j)}(y)</math>
thus in
:<math>y^i (1-y)^i M^{(i)}(y) = y^i (1-y)^i \mu^{(i)}(y) + y^i (1-y)^i s^{(i)} \ell(y) + y^i (1-y)^i \sum_{j=0}^{i-1} s^{(j)}(y) \ell^{(i-j)}(y)</math>

the first and third term are polynomials, and the middle term is a polynomial times a logit. Hence, this can always be written in the form:

:<math>y^i (1-y)^i M^{(i)}(y) = poly_1(y) + poly_2(y) \ell(y)</math>

where
:<math>poly_2(y) =  y^i (1-y)^i s^{(i)}</math>
:<math>poly_1(y) =  y^i (1-y)^i \mu^{(i)}(y) + T(y)</math>
:<math>T(y) = y^i (1-y)^i \sum_{j=0}^{i-1} s^{(j)}(y) \ell^{(i-j)}(y)</math>

When <math>i\ge\lfloor (k+1)/2 \rfloor</math>, then <math>poly_2(y)=0</math> so that <math>y^i (1-y)^i M^{(i)}(y)</math> is a simple polynomial.

=== Change of variable ===
<math>u = y-0.5</math>

<math>T(u) = T(y-0.5)</math>

== i<sup>th</sup> derivative finding algorithm ==

The algorithm uses arithmetic operations on a polynomial-coefficients data type. Here <math>T(u)</math>, <math>s^{(j)}</math>, <math>R(u)</math> and <math>N(u)</math> are all polynomials.  The order of the derivative, <math>i</math> is given.

#For each <math>j=0 \mbox{ to } i</math>
#* Compute <math>C_j = \left( \begin{array}{c}i  \\ j\end{array}\right) ( i-j-1)!</math>
#* Compute <math>s^{(j)}(u)</math>
#* Compute <math>R(u) = (0.5+u)^j (0.5-u)^j</math>
#* Compute <math>N(u) = (0.5+u)^{i-j} - (-1)^{i-j} ( 0.5-u)^{i-j}</math>
#* Compute <math>T_j(u) = C_j \cdot s^{(j)}(u) \cdot R(u) \cdot N(u)</math>
# <math>T(u) = \sum_{j=0}^i T_j(u)</math>

== Zero finding ==

To find the zeros of <math>M^{(i)}(y) </math>, first find the zeros of <math>M^{(i+1)}(y)</math>, which are the minima and maxima of <math>M^{(i)}(y) </math>. We should add to this list of extrema y=0 and y=1. Each adjacent extrema with opposite sign will bound one zero of <math>M^{(i)}(y) </math>. Thus, we have a finite set of intervals, each containing a single zero.

To find the zero, simply conduct a binary search within each interval that is known to contain the zero. Somewhat faster is a bounded Newton algorithm.