Relating to the article An approximation to the generalized hypergeometric function, I would like to calculated example of the Poisson distribution $Po(10)$. If you have the possibility, please see the article.
I am seeking advice on reproducing the example 5.1 of the poisson distribution $Po(10)$. According to table 2.1, (k,l,m,n)=(0,0,0,0). E.g. $N=1000$? $(a_0,b_0,c_0,d_0)=(1,1,1,1)$? And then $(\alpha_0,\beta_0,\gamma_0,\delta_0)=(0,0,0,0)$? But then is $A_2=B_2=0$.
I was trying to implement the formula (3.12) of Takeuchi's approximation approximation in Matlab (the results are different):
%% Takuchi (1984) approximation is formula (3.12)
lambda = 10; x=1;
trueVal=poisspdf(x,lambda);
N = 1000000; % N \to \infty
% eps=0.01;
% greek=eps/N; % greek=alpha=beta..
greek = 0;
theta=lambda; theta0=lambda/N; mu=N/lambda; % as stated in example 5.1
z=(x-N*mu)/sqrt(N);
B1 = 1/(mu)-1/(greek-mu)-1/(greek+mu)+1/(greek-mu);
sigma=1/sqrt(B1);
% sigma=mu;
B2 = 1/(greek+mu)^2-1/(greek-mu)^2-1/(greek+mu)^2+1/(greek-mu)^2;
B3 = 1/(greek+mu)^3-1/(greek-mu)^3-1/(greek+mu)^3+1/(greek-mu)^3;
A1 = 1/(greek+mu) +1/(greek-mu) -1/(greek+mu) -1/(greek-mu);
A2 = 1/(greek+mu)^2+1/(greek-mu)^2-1/(greek+mu)^2-1/(greek-mu)^2;
p = 1./(sqrt(2*pi*N)*sigma)*exp(-z^2/(2*sigma^2))*...
(...
1+1/(6*sqrt(N))*(A2*z^3-3*A1*z)-1/(24*N)*(2*B3*z^4-6*B2*z^2)+ ...
1/(72*N)*(A2*z^3-3*A1*z)^2 ...
);
relErr=(p-trueVal)/trueVal;
and the Engeworth approximation in formula (4.1)
%% Edgeworth approximation
% Example 5.1: Poisson Po(10)
lambda = 10; x=10 ; p=poisspdf(x,lambda);
sigma = lambda; mu = lambda;
n = lambda; % according to the article
y = x;
% z = ( y-n*mu+0.5 )/(sqrt(n)*sigma);
z = ( y-mu )/(sigma);
beta3 = mu/sigma^3; beta4 = mu/sigma^4;
pApprox = 1/( sqrt(2*pi*n*sigma) )*exp( -z^2/2 )*( 1+beta3/(6*sqrt(n))*( ...
z^3-3*z ) + beta4/(24*n)*( z^4-6*z^2+3 ) + beta3^2/(72*n)*( ...
z^6-15*z^4 ) );
relError = (p-pApprox)/p;
I will be happy for any suggestions!
Please see the formulas below (from the discussed article):
The generalizes hypergeometric distribution has the form GHG(k, l, m, n; $\vec{a}$, $\vec{b}$, $\vec{c}$, $\vec{d}$, $\theta$)
$$P\{X=x\}=c\cdot \frac{\prod_{j=1}^m c_j(x)\prod_{j=1}^n d_j^-[x]}{x!\prod_{j=1}^ka_j(x)\prod_{j=1}^l b_j^-[x]}\theta^x$$
For the Poisson case (k,l,m,n)=(0,0,0,0). Letting $\theta=\lambda$ and $\theta_0=\lambda/N$, we conclude $\mu=N/\lambda$
Using Takeuchi's approximation in the article, we should use $a_j=\alpha_jN+1$, $b_j=\beta_jN+1$, $c_j=\gamma_jN+1$, $d_j=\delta_jN+1$, $\alpha_j>0$, $\beta_j>0$, $\gamma_j>0$, for $\alpha_0=0$, we have $a_0=0$ (sure, whether also for the other variables).
We have
$$\begin{equation}\begin{split} A_i&=\sum_{j=0}^k\frac{1}{(\alpha_j+\mu)^i}+\sum_{j=1}^l\frac{1}{(\beta_j-\mu)^i}-\sum_{j=1}^m\frac{1}{(\gamma_j+\mu)^i}-\sum_{j=1}^n\frac{1}{(\delta_j-\mu)^i}\\ B_i&=\sum_{j=0}^k\frac{1}{(\alpha_j+\mu)^i}-\sum_{j=1}^l\frac{1}{(\beta_j-\mu)^i}-\sum_{j=1}^m\frac{1}{(\gamma_j+\mu)^i}+\sum_{j=1}^n\frac{1}{(\delta_j-\mu)^i} \end{split}\end{equation}$$ which we input into \begin{equation}\label{eq:approxHida} \begin{split} p_X(x)&=\frac{1}{\sqrt{2\pi N\sigma}}e^{-\frac{z^2}{2\sigma^2}}\Big[1+(A_2z^3-3A_1z)-\frac{1}{24N}( 2B_2z^4-6B_2z^2)\\ &+\frac{1}{72N}(A_2z^3-3A_1z)^2+o\left(\frac{1}{N}\right)\Big], \end{split} \end{equation} where $z=(x-N\mu)/\sqrt{N}$.
For the Engeworth approximation:
$$P\{Y_n=y\}=\frac{1}{\sqrt{2\pi n}\sigma}e^{-z^2/2}\Big( 1+\frac{\beta_3}{6\sqrt{n}}H_3(z)+\frac{\beta_4}{24\sqrt{n}}H_3(z)+\frac{\beta_3^2}{72\sqrt{n}}H_6(z)\Big)$$, where $z=(y-n\mu+0.5)/(\sqrt{n}\sigma)$ $H_i$'s the Hermite polynomials: $H_3(z)=z^3-3z$, $H_4(z)=z^4-6z^2+3$, $H_6(z)=z^6-15z^4+45z^2-15$. $\beta_i=\kappa_i/\sigma^i$. $\kappa_i$ is the i-th cummulant. (in the poisson example are the $\beta_i=\frac{\lambda}{\lambda^i}$ and $\lambda=n$).
The results should be for $x=1$, the true value 0.00045399 and the relative error for Takeuchi -0.112287 and Edgeworth. -0.163343. The results should be for $x=2$, the true value 0.00227 and the relative error for Takeuchi 0.0604028 and Edgeworth -0.163343.