Julia - double integral

Question

I have two functions written in Julia that I want to integrate and both should give a value very close to 1. They both produce the same type of variables but one integral works and the other doesn't. Can someone help me?

The one that works

using Cubature
function joedensity(u,v,alpha)
    w = 1-u
    z = 1-v
    (w^alpha+z^alpha-(w*z)^alpha)^(1/alpha-2)*(w*z)^(alpha-1)*(alpha-1+w^alpha+z^alpha-(w*z)^alpha)
end

u = [0.1, 0.5, 0.3, 0.7, 0.9]
v = [0.01, 0.2, 0.8, 0.3, 0.99]
alpha = 3
joedensity.(u,v,alpha)
#= Results
5-element Vector{Float64}:
 2.407462338226518
 1.0414427450537544
 0.2669954330753848
 0.5695056921157302
 0.1997673764704445
=#
hcubature(x -> joedensity.(x[1],x[2],alpha),[0,0],[1,1],abstol=0)
#= Results
(1.0000000000808127, 9.994730188502498e-9)
=#

The one that doesn't

this function calls other functions that I defined, but it is working exactly as I want

function quanty(q,delta)
    dummyd = y -> (1-q)-(delta/(2*delta-1))*exp(-y/delta)+((1-delta)/(2*delta-1))*exp(-y/(1-delta))
    dummy05 = y -> (1-q)-exp(-2*y)*(2*y+1)
    (abs(delta-0.5)>1e-05) ? find_zero(dummyd,(0,100)) : find_zero(dummy05,(0,100))
end

function jgaus(s,thetaW)
    z = map(q -> quantile.(Normal(0,1), 1 .-exp.(-q)), s)
    if any(z==Inf)
        z = -map(q -> quantile.(Normal(0,1), exp.(-q)), s)
    end
    inside = exp.(-(1/(2*(1-thetaW^2)))*(z[1].^2 - 2*thetaW*z[1].*z[2] + z[2].^2))
    (1/(2*pi*sqrt(1-thetaW^2)))*inside .* prod(map(q -> exp.(-q), s) ./ map(q -> pdf.(Normal(0,1), q), z))
end

function jointdist(y,delta,thetaW)
    function int(v)
        s = map(q -> (q .- delta*v) ./ (1-delta), y)
        jgaus(s, thetaW) .* exp.(-v)
    end
    integ, err = quadgk(int,0,minimum(y)./delta,atol=0)
    integ*(1-delta)^(-2)
end

marg = (y,delta) -> (abs(delta-0.5)>1e-05) ? ((1/(2*delta-1))*exp(y)^(-1/delta-1)-(1/(2*delta-1))*exp(y)^(-1/(1-delta)-1))*exp(y) : ((exp(y)^(-3)*(4*y)))*exp(y)

using Cubature
function cop2019(u,v,delta,thetaW)
    y = zeros(2,length(u))
    y[1,:] .= quanty.(u,delta)
    y[2,:] .= quanty.(v,delta)
    map(x -> jointdist(x,delta,thetaW), eachcol(y)) ./ (marg.(y[1,:],delta) .* marg.(y[2,:],delta))
end 

u = [0.1, 0.5, 0.3, 0.7, 0.9]
v = [0.01, 0.2, 0.8, 0.3, 0.99]
delta = 0.2
thetaW = 0.4;
cop2019(u,v,delta,thetaW)
#= Results
5-element Vector{Float64}:
 2.044469250947409
 0.9381167412987401
 0.7335531873911258
 0.840984518753334
 2.3684851388756365
=#
hcubature(x -> cop2019(x[1],x[2],delta,thetaW),[0,0],[1,1],abstol=0)
#= Part of the error
MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, ::Base.TwicePrecision) where T<:Number at ~/julia-1.7.2/share/julia/base/twiceprecision.jl:262
  convert(::Type{T}, ::AbstractChar) where T<:Number at ~/julia-1.7.2/share/julia/base/char.jl:185
  convert(::Type{T}, ::CartesianIndex{1}) where T<:Number at ~/julia-1.7.2/share/julia/base/multidimensional.jl:136
  ...
=#

Answer 1

The problem becomes apparent when you call the functions the way hcubature calls them, with individual values:

julia> cf = x -> cop2019(x[1],x[2],delta,thetaW);

julia> jf = x -> joedensity.(x[1],x[2],alpha);

julia> xin = [0.1, 0.2];

julia> cf(xin)
1-element Vector{Float64}:
 1.5478771875651165

julia> jf(xin)
1.883131054500218

The basic difference is that joedensity is called with a broadcast, and when such a call is made with only scalar values, its return value is also a scalar. cop2019 however always returns a Vector value because that's how the function has been defined.

A temporary fix for this is to change its last line to only(map(x -> jointdist(x,delta,thetaW), eachcol(y)) ./ (marg.(y[1,:],delta) .* marg.(y[2,:],delta))).

Trying that change, the original type error about convert goes away; there's a new DomainError about the integrand function, but I believe that error is more domain-specific (no pun intended) to your numerical problem.

ERROR: DomainError with 0.5015637133281883:
integrand produced NaN in the interval (0.5015637133281599, 0.5015637133282168)

However, a more correct fix would be to change the cop2019 function to work on individual values rather than whole vectors, and just like joedensity, have it return a single value.