ERROR: LoadError: DomainError with 0.5: integrand produced NaN in the interval (0.0, 1.0)

[aaabdulhaq]

Following is my code to solve two-dimensional integral in Julia.

using Statistics,Distributions,Distributed,StatsFuns,Roots,QuadGK
p=3;n0=3;n1=1;n2=4;d=0.25;
inARL=200;n=n1+n2;w1=chisqinvcdf(p,1-(n0-n1)/n2);inpa1=chisqcdf(p,w1);

function inpa(h2)
function fvv(v)
return(nchisqcdf(p,(n1/n2)v,(n*h2)/n2)*chisqpdf(p,v));end;
return(1/(1-inpa1-quadgk(fvv,w1,Inf,rtol=1e-10)[1])-inARL);end;

h=fzero(inpa,0.25,50);
outpa1=nchisqcdf(p,n1*d^2,w1);

function fz1(z)
function fv1(v)
lam=(n1/n2)*((z+sqrt(n1)*d)^2+v)+(2*n*sqrt(n1)*d*z)/(n2)+((n*d)^2)/(n2);
return(nchisqcdf(p,lam,(n*h)/n2)*chisqpdf(p-1,v));end;
return(normpdf(z)*quadgk(fv1,0,Inf,rtol=1e-10)[1]);end;
int1=quadgk(fz1,-Inf,-sqrt(w1)-sqrt(n1)*d,rtol=1e-10)[1];# problem

[stevengj]

Apparently your integrand returned a NaN? Have you checked that your integrands work in the intervals you specify?

(I should also note that for 2d integrals you are usually better off with a 2d integration routine like HCubature.jl than nesting 1d adaptive quadratures.)

[aaabdulhaq]

Thank you for your response. Yes, the integral converges on the ranges as I've computed the same result in Mathematica. I don't know how to implement 2d integration in HCubature as they require all integrals' ranges to be in (0,1). In my case, the integrals are dependent and indefinite.

[stevengj]

Yes, the integral converges on the ranges as I've computed the same result in Mathematica.

The error message indicates that your integrand function is returning NaN for some values of the argument. This may indicate a bug in your code.

I don't know how to implement 2d integration in HCubature as they require all integrals' ranges to be in (0,1). In my case, the integrals are dependent and indefinite.

You use a change of variables as described here: https://github.com/stevengj/cubature#infinite-intervals

(This is what QuadGK is doing internally.)

[aaabdulhaq]

I followed you but the problem is still there! Following is the simplest code with change of variable implemented.

using Statistics,Distributions,Distributed,StatsFuns,Roots,QuadGK
p=3;d=0.25;h=13.84283;n1=1;n2=4;n=n1+n2;

function fz1(z)
function fv1(v)
lam=(n1/n2)*(((z+sqrt(n1)*d)^2)+v/(1-v))+(2*n*sqrt(n1)*d*z)/(n2)+((n*d)^2)/(n2);
return(nchisqcdf(p,lam,(n*h)/n2)*chisqpdf(p-1,v/(1-v))*(1/((1-v)^2)));end;
return(quadgk(fv1,0,1,rtol=1e-10)[1]);end;

fz1(0) # okay
fz1(-1)   # problem
fz1(-2)   # problem
fz1(-3)   # okay

julia> fz1(-3)
0.9964662237851157

julia> fz1(-2)
ERROR: DomainError with 0.5:
integrand produced NaN in the interval (0.0, 1.0)
Stacktrace:
 [1] evalrule(::getfield(Main, Symbol("#fv1#3")){Int64}, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::typeof(LinearAlgebra.norm)) at C:\Users\AH\.julia\packages\QuadGK\fTRFu\src\QuadGK.jl:81
 [2] do_quadgk(::getfield(Main, Symbol("#fv1#3")){Int64}, ::Array{Float64,1}, ::Int64, ::Type{Float64}, ::Float64, ::Float64, ::Int64, ::typeof(LinearAlgebra.norm)) at C:\Users\AH\.julia\packages\QuadGK\fTRFu\src\QuadGK.jl:124
 [3] #quadgk at C:\Users\AH\.julia\packages\QuadGK\fTRFu\src\QuadGK.jl:170 [inlined]
 [4] #quadgk#15 at C:\Users\AH\.julia\packages\QuadGK\fTRFu\src\QuadGK.jl:244 [inlined]
 [5] #quadgk at .\none:0 [inlined]
 [6] fz1(::Int64) at d:\julia\int.jl:57
 [7] top-level scope at none:0

julia> fz1(-1)
ERROR: DomainError with 0.5:
integrand produced NaN in the interval (0.0, 1.0)
Stacktrace:
 [1] evalrule(::getfield(Main, Symbol("#fv1#3")){Int64}, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::typeof(LinearAlgebra.norm)) at C:\Users\AH\.julia\packages\QuadGK\fTRFu\src\QuadGK.jl:81
 [2] do_quadgk(::getfield(Main, Symbol("#fv1#3")){Int64}, ::Array{Float64,1}, ::Int64, ::Type{Float64}, ::Float64, ::Float64, ::Int64, ::typeof(LinearAlgebra.norm)) at C:\Users\AH\.julia\packages\QuadGK\fTRFu\src\QuadGK.jl:124
 [3] #quadgk at C:\Users\AH\.julia\packages\QuadGK\fTRFu\src\QuadGK.jl:170 [inlined]
 [4] #quadgk#15 at C:\Users\AH\.julia\packages\QuadGK\fTRFu\src\QuadGK.jl:244 [inlined]
 [5] #quadgk at .\none:0 [inlined]

julia> fz1(-3)
0.9964662237851157

julia> fz1(-0)
0.9964662237851157

[aaabdulhaq]

Now I compare these results with Mathematica code:

In[4]:= p = 3; n1 = 1; n2 = 4; n = n1 + n2; d = 0.25; h = 13.84;

In[5]:= INT1[
   z_?NumericQ] := (Return[
     NIntegrate[
      CDF[NoncentralChiSquareDistribution[p, 
         n1/n2*((z + Sqrt[n1]*d)^2 + v) + (2.*n*Sqrt[n1]*d*z)/
          n2 + (n*d)^2/n2], (n*h)/n2]*
       PDF[ChiSquareDistribution[p - 1], v][[1]][[1]][[1]], {v, 
       0., \[Infinity]}]];);

In[6]:= {INT1[-3], INT1[-2], INT1[-1], INT1[0]}

Out[6]= {0.996462, 0.998166, 0.998166, 0.996462}

[aaabdulhaq]

The results are of course okay for the arguments -3 and 0; but, for others, Julia refuses to answer and reports NaN. But, Mathematica reports the answer.

[stevengj]

The problem is that your integrand is returning NaN for e.g. z=-1 and v=0.1, exactly as it says:

using StatsFuns
p=3;d=0.25;h=13.84283;n1=1;n2=4;n=n1+n2;
function fv(z, v)
    lam=(n1/n2)*(((z+sqrt(n1)*d)^2)+v/(1-v))+(2*n*sqrt(n1)*d*z)/(n2)+((n*d)^2)/(n2)
    return(nchisqcdf(p,lam,(n*h)/n2)*chisqpdf(p-1,v/(1-v))*(1/((1-v)^2)))
end

gives

julia> fv(-1, 0.1)
NaN

in particular, the nchisqcdf function is returning NaN. With the help of a couple of @show statements, I can see that you are evaluating:

julia> nchisqcdf(3, -0.06597222222222221, 17.303537499999997)
NaN

It looks like nchisqcdf doesn't like any negative value for the second argument. If you think this is wrong, you could file an issue with the StatsFuns.jl package.

However, in JuliaStats/StatsFuns.jl#19 it was stated that you are not supposed to call nchisqcdf directly (it is undocumented), and are supposed to use Distributions.jl instead.