Search code examples
matlabdoubleintegralequation-solving

Solving a complex double-integral equation for a third variable

I'm trying to solve the following equation for the variable h: enter image description here

  • F is the normal cumulative distribution function
  • f is the probability density function of chi-square distribution.

I want to solve it numerically using Matlab.

First I tried to solve the problem with a simpler version of function F and f.

Then, I tried to solve the problem as below:

n = 3;    % n0 in the equation
k = 2;    
P = 0.99; % P* in the equation
F = @(x,y,h) normcdf(h./sqrt((n-1)*(1./x+1./y)));
g1 = @(y,h)integral(@(x) F(x,y,h).*chi2pdf(x,n),0,Inf, 'ArrayValued', true);
g2 = @(h) integral(@(y) (g1(y,h).^(k-1)).*chi2pdf(y,n),0,Inf)-P;
bsolve = fzero(g2,0)

I obtained an answer of 5.1496. The problem is that this answer seems wrong.

There is a paper which provides a table of h values obtained by solving the above equation. This paper was published in 1984, when the computer power was not good enough to solve the equation quickly. They solved it with various values:

  • n0 = 3:20, 30:10:50
  • k = 2:10
  • P* = 0.9, 0.95 and 0.99

They provide the value of h in each case.

Now I'm trying to solve this equation and get the value of h with different values of n0, k and P* (for example n0=25, k=12 and P*=0.975), where the paper does not provide the value of h.

The Matlab code above gives me different values of h than the paper. For example:

  • n0=3, k=2 and P*=0.99: my code gives h = 5.1496, the paper gives h = 10.276.
  • n0=10, k=4 and P*=0.9: my code gives h = 2.7262, the paper gives h = 2.912.

The author says

The integrals were evaluated using 32 point Gauss-Laguerre numerical quadrature. This was accomplished by evaluating the inner integral at the appropriate 32 values of y which were available from the IBM subroutine QA32. The inner integral was also evaluated using 32 point Gauss-Laguerre numerical quadrature. Each of the 32 values of the inner integral was raised to the k-1 power and multiplied by the appropriate constant.

Matlab seems to use the same method of quadrature:

q = integral(fun,xmin,xmax) numerically integrates function fun from xmin to xmax using global adaptive quadrature and default error tolerances.

Does any one have any idea which results are correct? Either I have made some mistake(s) in the code, or the paper could be wrong - possibly because the author used only 32 values in estimating the integrals?

Solution

  • This does work, but the chi-squared distribution has (n-1) degrees of freedom, not n as in the code posted. If that's fixed then the Rinott's constant values agree with literature. Or at least, the literature that isn't behind a paywall. Can't speak to the 1984 Wilcox.