kappa34.mws

kappa34.mws  030128

Determine cumulants for 1) Poisson distribution,  2) Binomial distribution,  3) Log-normal distribution

>    restart;

1) Begin with Poisson distribution. Determine the first 4 moments about zero.

>    sum(lambda^k*exp(-lambda)/k!,k=0..infinity);

>    alpha1:=sum(k*lambda^k*exp(-lambda)/k!,k=0..infinity);

>    alpha2:=sum(k^2*lambda^k*exp(-lambda)/k!,k=0..infinity);

>    alpha3:=sum(k^3*lambda^k*exp(-lambda)/k!,k=0..infinity);

>    alpha4:=sum(k^4*lambda^k*exp(-lambda)/k!,k=0..infinity);

1

alpha1 := lambda

alpha2 := lambda*(lambda+1)

alpha3 := lambda*(1+3*lambda+lambda^2)

alpha4 := lambda*(1+7*lambda+6*lambda^2+lambda^3)

The cumulants are determined from the moments:

>    kappa1:=alpha1;
kappa2:=expand(alpha2-alpha1^2);
kappa3:=expand(alpha3-3*alpha2*alpha1+2*alpha1^3);
kappa4:=expand(alpha4-4*alpha3*alpha1-3*alpha2^2+12*alpha2*alpha1^2-6*alpha1^4);

kappa1 := lambda

kappa2 := lambda

kappa3 := lambda

kappa4 := lambda

2) Now turn to a binomially distributed r.v.:

>    restart;

The first four moments are determined as functions of the parameters N and p:

>    simplify(sum(binomial(N,k)*p^k*(1-p)^(N-k),k=0..N),symbolic);

>    alpha1:=simplify(sum(k*binomial(N,k)*p^k*(1-p)^(N-k),k=0..N),symbolic);

>    alpha2:=map(simplify,sum(k^2*binomial(N,k)*p^k*(1-p)^(N-k),k=0..N),symbolic);

>    alpha3:=map(simplify,sum(k^3*binomial(N,k)*p^k*(1-p)^(N-k),k=0..N),symbolic);

>    alpha4:=map(simplify,sum(k^4*binomial(N,k)*p^k*(1-p)^(N-k),k=0..N),symbolic);

1

alpha1 := N*p

alpha2 := N*p+p^2*N*(N-1)

alpha3 := N*p+3*p^2*N*(N-1)+p^3*N*(N-1)*(N-2)

alpha4 := N*p+7*p^2*N*(N-1)+6*p^3*N*(N-1)*(N-2)+p^4*N*(N-1)*(N-2)*(N-3)

The cumulants are determined from the moments:

>    kappa1:=alpha1;

>    kappa2:=expand(alpha2-alpha1^2);

>    kappa3:=expand(alpha3-3*alpha2*alpha1+2*alpha1^3);

>    kappa4:=expand(alpha4-4*alpha3*alpha1-3*alpha2^2+12*alpha2*alpha1^2-6*alpha1^4);

kappa1 := N*p

kappa2 := N*p-N*p^2

kappa3 := N*p-3*N*p^2+2*N*p^3

kappa4 := -7*N*p^2+N*p+12*N*p^3-6*p^4*N

I wish to express kappa3 and kappa4 in terms of kappa1 and kappa2. Solve the first two equations above for N and p in terms of kappa1 and kappa2:  Get kappa1*(1-p)=kappa2; hence p=1-kappa2/kappa1.  

 

>    simplify(subs(N=kap1/(1-kap2/kap1),p=1-kap2/kap1,kappa3));

-1/kap1*kap2*(kap1-2*kap2)

>    simplify(subs(N=kap1/(1-kap2/kap1),p=1-kap2/kap1,kappa4));

1/kap1^2*kap2*(kap1^2-6*kap1*kap2+6*kap2^2)

3) Log-normal distribution: Y=ln(X/x0) is N(mu,sigma). I have X=x0*exp(Y). Hence

   EX^n  =   int(x0^n*exp(n*y)*phi((y-mu)/sigma)/sigma,y = -infinity .. infinity)  , where phi  is the normal density. By completing the square and integrating, I get

  alpha_n  = x0^n*exp(n*mu+n^2*sigma^2/2)   .   

>    restart;

The moments are determined from mu and sigma as folloes:

>    alpha:=n->x0^n*exp(n*mu+n^2*sigma^2/2);

alpha := proc (n) options operator, arrow; x0^n*exp(n*mu+1/2*n^2*sigma^2) end proc

The first four cumulants:

>    kappa1:=alpha(1);

>    kappa2:=simplify(alpha(2)-alpha(1)^2);

>    kappa3:=simplify(alpha(3)-3*alpha(1)*alpha(2)+2*alpha(1)^3);

>    kappa4:=simplify(alpha(4)-3*alpha(2)^2-4*alpha(1)*alpha(3)+12*alpha(1)^2*alpha(2)-6*alpha(1)^4);

kappa1 := x0*exp(mu+1/2*sigma^2)

kappa2 := -x0^2*(-exp(2*mu+2*sigma^2)+exp(2*mu+sigma^2))

kappa3 := x0^3*(exp(3*mu+9/2*sigma^2)-3*exp(3*mu+5/2*sigma^2)+2*exp(3*mu+3/2*sigma^2))

kappa4 := x0^4*(exp(4*mu+8*sigma^2)-3*exp(4*mu+4*sigma^2)-4*exp(4*mu+5*sigma^2)+12*exp(4*mu+3*sigma^2)-6*exp(4*mu+2*sigma^2))

Note that kappa_k contains the factor kappa_1^k. I simplify the expressions for the cumulants:

>    expand(simplify(kappa2/kappa1^2,power));

>    simplify(expand(simplify(kappa3/kappa1^3,power)),power);

>    simplify(expand(simplify(kappa4/kappa1^4,power)),power);

exp(sigma^2)-1

exp(3*sigma^2)-3*exp(sigma^2)+2

exp(6*sigma^2)-3*exp(2*sigma^2)-4*exp(3*sigma^2)+12*exp(sigma^2)-6

I wish to express the cumulants kappa3 and kappa4 in terms of kappa1 and kappa2.

The first of these expressions is solved for exp(sigma^2) in terms of kappa2 and kappa1. I get exp(sigma^2)=1+kappa2/kappa1^2. Use the notation kap_k for kappa_k, and

>    esig:=1+kap2/kap1^2;

esig := 1+kap2/kap1^2

>    kap3:=expand(kap1^3*(esig^3-3*esig+2));

kap3 := 3/kap1*kap2^2+1/kap1^3*kap2^3

>    kap4:=expand(kap1^4*(esig^6-4*esig^3-3*esig^2+12*esig-6));

kap4 := 16/kap1^2*kap2^3+15/kap1^4*kap2^4+6/kap1^6*kap2^5+1/kap1^8*kap2^6