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

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

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

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

1 for i = 0 to n loop
2 wi ← 1 // wi = φ0(t; x1, x2, …, xi) = 1
3 end loop
4 for k = 1 toloop
5 w0w0t / k // w0 = φk(t; ~) = tk / k!
6 for i = 1 to n loop // O(n) loop
7 wiwi − 1 + xiwi // wi = φk(t; x1, x2, …, xi)
8 end loop
9 output wn // wn = φk(t; x1, x2, …, xn)
10 end 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 (gDμj)(t), and the coefficients Ai, k change as shown below:

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