朱莉娅 - 二重积分

Question

我有两个用 Julia 编写的函数，我想要对其进行积分，并且两个函数都应该给出非常接近 1 的值。它们都生成相同类型的变量，但一个积分有效，而另一个则无效。有人可以帮助我吗？

有效的

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)
=#

无效的

这个函数调用我定义的其他函数，但它完全按照我想要的方式工作

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

当您以 hcubature 调用函数的方式（使用单独的值）调用函数时，问题就变得很明显：

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

基本区别是 joedensis 是通过广播调用的，并且当这样的调用时仅由标量值组成，其返回值也是标量。然而，cop2019 始终返回一个 Vector 值，因为这就是函数的定义方式。

对此的临时修复是将其最后一行更改为 only(map(x -> jointdist(x,delta,thetaW),eachcol(y)) ./ (marg.(y[1,:], delta) .* marg.(y[2,:],delta)))。

尝试进行更改后，有关 convert 的原始类型错误消失了；关于被积函数有一个新的 DomainError ，但我相信该错误对于您的数值问题来说更特定于域（没有双关语）。

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

然而，更正确的修复方法是更改​​ cop2019 函数以处理单个值而不是整个向量，并且就像 joedensis 一样，让它返回单个值。