求解第三个变量的复数二重积分方程
Solving a complex double-integral equation for a third variable
我正在尝试为变量 h 求解以下等式：
- F是正态累积分布函数
- f是卡方分布的概率密度函数
我想用Matlab数值求解。
首先我尝试。
然后，我尝试解决问题如下：
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)
我得到了5.1496的答案。问题是这个答案似乎是错误的。
有一篇论文提供了table解上述方程得到的h值。这篇论文发表于 1984 年，当时计算机的能力还不足以快速求解方程。他们用各种值解决了它：
- n0 =
3:20, 30:10:50
- k =
2:10
- P* =
0.9、0.95 和 0.99
它们提供了每种情况下的 h 值。
现在我正在尝试求解这个方程并得到具有不同 n0、k 和 P* 值的 h 值（例如 n0=25、k=12 和 P*=0.975），其中论文不提供 h 的值。
上面的 Matlab 代码给出了与论文不同的 h 值。
例如：
- n0=3, k=2 and P*=0.99: 我的代码给出h = 5.1496，论文给出h = 10.276。
- n0=10, k=4 and P*=0.9: 我的代码给出h = 2.7262，论文给出h = 2.912.
作者说
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好像用的求积法一样：
q = integral(fun,xmin,xmax) numerically integrates function fun from xmin to xmax using global adaptive quadrature and default error tolerances.
有人知道哪些结果是正确的吗？要么我在代码中犯了一些错误，要么这篇论文可能是错误的——可能是因为作者在估计积分时只使用了 32 个值？
这确实有效，但卡方分布具有 (n-1) 个自由度，而不是发布的代码中的 n。如果这是固定的，那么 Rinott 的常数值与文献一致。或者至少，没有付费专区的文献。无法与 1984 年的 Wilcox 对话。
我正在尝试为变量 h 求解以下等式：
- F是正态累积分布函数
- f是卡方分布的概率密度函数
我想用Matlab数值求解。
首先我尝试
然后，我尝试解决问题如下：
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)
我得到了5.1496的答案。问题是这个答案似乎是错误的。
有一篇论文提供了table解上述方程得到的h值。这篇论文发表于 1984 年，当时计算机的能力还不足以快速求解方程。他们用各种值解决了它：
- n0 =
3:20,30:10:50 - k =
2:10 - P* =
0.9、0.95和0.99
它们提供了每种情况下的 h 值。
现在我正在尝试求解这个方程并得到具有不同 n0、k 和 P* 值的 h 值（例如 n0=25、k=12 和 P*=0.975），其中论文不提供 h 的值。
上面的 Matlab 代码给出了与论文不同的 h 值。 例如：
- n0=3, k=2 and P*=0.99: 我的代码给出h = 5.1496，论文给出h = 10.276。
- n0=10, k=4 and P*=0.9: 我的代码给出h = 2.7262，论文给出h = 2.912.
作者说
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好像用的求积法一样：
q = integral(fun,xmin,xmax) numerically integrates function fun from xmin to xmax using global adaptive quadrature and default error tolerances.
有人知道哪些结果是正确的吗？要么我在代码中犯了一些错误，要么这篇论文可能是错误的——可能是因为作者在估计积分时只使用了 32 个值？
这确实有效，但卡方分布具有 (n-1) 个自由度，而不是发布的代码中的 n。如果这是固定的，那么 Rinott 的常数值与文献一致。或者至少，没有付费专区的文献。无法与 1984 年的 Wilcox 对话。