Algorithms

Definitions
chspHR.png (1)
chsp-hkmxHR.png (2)
decay-phi-defHR.png
(3)

1. Incremental calculation of hk(x1, x2, …, xn), for increasing k, with O(n) increment step

1for i = 1 to n loop
2wi ← 1 // wi = h0(x1, x2, …, xi) = 1
3end loop
4for k = 1 toloop
5 w1x1w1 // w1 = hk(x1) = x1k
6for i = 2 to n loop
7wiwi − 1 + xiwi // wi = hk(x1, x2, …, xi)
8end loop
9output wn // wn = hk(x1, x2, …, xn)
10end loop

2. Incremental calculation of hk(m1x1, m2x2, …, mnxn) — which may be useful when ∑ mi n

1y ← 1 // y = h0(m1x1, m2x2, …, mnxn) = 1
2for i = 1 to n loop  
3wi ← 0  
4end loop  
5for k = 1 toloop  
6s ← 0  
7for i = 1 to n loop // O(n) loop
8wixi (wi + miy)  
9ss + wi  
10end loop  
11ys / k // y = hk(m1x1, m2x2, …, mnxn)
12output y  
13end loop  

The same algorithm may be useful when the mi are not integers (e.g., when calculating a convolution of gamma distributions with varying scale parameters). In this case the algorithm outputs for each k the value

chsp-yHR.png

3. Incremental calculation of φk(t; x1, x2, …, xn)

1for i = 0 to n loop
2wi ← 1 // wi = φ0(t; x1, x2, …, xi) = 1
3end loop
4for k = 1 toloop
5w0w0t / k // w0 = φk(t; ~) = tk / k!
6for i = 1 to n loop // O(n) loop
7wiwi − 1 + xiwi // wi = φk(t; x1, x2, …, xi)
8end loop
9output wn // wn = φk(t; x1, x2, …, xn)
10end loop

4. If G(s) is the rational function

chsp-GsHR.png

then the inverse Laplace transform of G(s) is the function

chsp-gt1HR.png

where Ai, k is given by:

decay-Aik-2HR.png

If one modifies G(s) by incrementing mj for some 1 ≤ jn, then g(t) is replaced by (g Dμj)(t), and the coefficients Ai, k change as shown below:

1s ← 0
2for i = 1 to n loop
3if ij then
4Ai, 0Ai, 0 / (μjμi)
5for k = 1 to mi − 1 loop
6Ai, k ← (Ai, kAi, k − 1) / (μjμi)
7end loop
8ss + Ai, mi − 1
9end if
10end loop
11Aj, mj ← −s // new coefficient
12mjmj + 1

If j = n + 1 so that we are including a new factor (s + μn + 1) in the denominator of G(s), insert the line mj ← 0 before line 1 and append the line nn + 1 after line 12.

If the value of g(t) is needed only for a specific value of t, and if it is known that one mj will be incremented repeatedly and no other mi will increase, then a different version of the algorithm can be used to recalculate Ai, k only for ij, without accumulating the sum s. The sum of all the Aj, k terms can be recalculated at each iteration using the φ function above without calculating the individual coefficients Aj, k.

Note: Given any function g(t) of the form

chsp-gt2HR.png

with arbitrary coefficients Ai, k, if we wish to calculate (g Dμj)(t), the same algorithm for updating the coefficients works. However, the shortcut method of evaluating g(t) using the φ function is not guaranteed to work for arbitrary Ai, k.