Sunday, April 17, 2016

HP Prime: Gompertz Model

HP Prime:  Gompertz Model

The program GOMPERTZ fits the data {y_1, y_2, y_3, … , y_m } to the curve:

y = a*b^c^x

where x = {0, 1, 2, 3, 4, … , m-1}.  The program works best if the dimension (size) of the data set is divisible by 3.  GOMPERTZ sorts the values into ascending order before calculation.

The results returned is a 3 x 1 matrix of the parameters: a, b, and c.    

GOMPERTZ estimates a, b, and c and runs through one iteration to retrieve a better estimate, and may be modified by include additional iterations should the user desire. 

Process (in a nutshell):

1.  Let y = list of data to fit.  Sort y to list all elements in ascending order.

2.  Let m = size(y).  Divide y into three equal-sized sub-lists.   For each sub-list, take the natural logarithm of each element and add.  Hence:
n = m/3
s1 = Σ(ln y_k) for k=1 to n
s2 = Σ(ln y_k) for k=n+1 to 2n
s3 = Σ(ln y_k) for k=2n+1 to m

3.  Calculate the first estimate of a, b, and c:
c = ((s3 – s2)/(s2 –s1))^(1/n)
a = e^( 1/n * (s1 + (s2 – s1)/(1 – c^n))
b = e^( (s2 – s1)*(c – 1)/(1 – c^n)^2)

4.  Create the matrix v, where v = [[a],[b],[c]]. 

5.  Develop the matrix D and Y, where:

Each row of D consists of
[ b^c^I, b^c^I * a/b * c^I,  b^c^I * a/b * c^I * ln(b)* I]

And each row of Y consists of
[ y(I+1) – a*b^c^I ]

Where I = 0 to m-1

6.  Calculate v’ by setting v’ = (D^T*D)^-1 * D^T * Y. 

7.  Add v’ to v.  v1 = v + v’.  The matrix v1 is the new estimate.

The equation is estimated to be:  y ≈ a/100 * b^c^x   

HP Prime Program: GOMPERTZ

EXPORT GOMPERTZ(yl)
BEGIN
LOCAL y,m,n;
// 2nd STEP
y:=SORT(yl);
m:=SIZE(y);
n:=m/3;
LOCAL y1:=SUB(y,1,n);
LOCAL y2:=SUB(y,n+1,2*n);
LOCAL y3:=SUB(y,2*n+1,m);
LOCAL s1:=ΣLIST(LN(y1));
LOCAL s2:=ΣLIST(LN(y2));
LOCAL s3:=ΣLIST(LN(y3));
LOCAL c:=((s3-s2)/(s2-s1))^(1/n);
LOCAL b:=e^((s2-s1)*(c-1)/
(1-c^n)^2);
LOCAL a:=e^(1/n*(s1+(s2-s1)/
(1-c^n)));
LOCAL v:=[[a],[b],[c]];

// Making dmat,ymat
LOCAL dmat,ymat,I,k1,k2,k3,k4;
dmat:=[[0,0,0]];
ymat:=[[0]];
FOR I FROM 0 TO m-1 DO
k1:=b^c^I;
k2:=k1*a/b*c^I;
k3:=k2*LN(b)*I;
k4:=y(I+1)-a*(b^c^I);
dmat:=ADDROW(dmat,[k1,k2,k3],
I+1);
ymat:=ADDROW(ymat,[k4],I+1);
END;
dmat:=DELROW(dmat,I+1);
ymat:=DELROW(ymat,I+1);

// new param
LOCAL vt:=(TRN(dmat)*dmat)*
TRN(dmat)*ymat;
LOCAL v1:=v+vt;
RETURN v1;

END;


Example

Data: {58, 66, 72.5, 78, 82, 85}

Result:

[[ 94.2216370902 ] , [ 0.615221033606 ] , [ 0.732101473203 ]]

We can infer that curve is estimated to be:

y ≈ 94.2216370902 * 0.615221033606^0.732101473203^x


Source

ReliaWiki.  “Gompertz Models”  http://reliawiki.org/index.php/Gompertz_Models   Retrieved April 12, 2016



This blog is property of Edward Shore, 2016