Algorithms

Definitions
chsp.png (1)
chsp-hkmx.png (2)
decay-phi-def.png
(3)

1. Incremental calcula­tion of hk(x1, x2, …, xn), for increas­ing k, with O(n) incre­ment 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 calcula­tion of hk(m1x1, m2x2, …, mnxn) — which may be use­ful 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 algo­rithm may be use­ful when the mi are not inte­gers (e.g., when cal­cu­lating a con­volu­tion of gamma dis­tribu­tions with vary­ing scale param­eters). In this case the algo­rithm out­puts for each k the value

chsp-y.png

3. Incremental calcula­tion 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-Gs.png

then the inverse Laplace trans­form of G(s) is the function

chsp-gt1.png

where Ai, k is given by:

decay-Aik-2.png

If one modifies G(s) by incre­ment­ing mj for some 1 ≤ jn, then g(t) is replaced by (g Dμj)(t), and the co­effi­cients 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 fac­tor (s + μn + 1) in the denom­inator of G(s), in­sert 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 incre­mented repeatedly and no other mi will in­crease, then a dif­fer­ent ver­sion of the algo­rithm can be used to recal­cu­late Ai, k only for ij, with­out accumu­lating the sum s. The sum of all the Aj, k terms can be recal­cu­lated at each itera­tion using the φ function above with­out cal­cu­lating the in­di­vid­ual co­effi­cients Aj, k.

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

chsp-gt2.png

with arbitrary coeffi­cients Ai, k, if we wish to cal­cu­late (g Dμj)(t), the same algo­rithm for up­dating the co­effi­cients works. However, the short­cut method of evalu­ating g(t) using the φ func­tion is not guar­an­teed to work for arbi­trary Ai, k.