Hcubature - New to Julia

New to Julia
1

Hi all! I’m very new to Julia so this may be quite simple I don’t know.
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

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
  ...
=#
2

It expects your function to return a scalar, but it is returning a vector.

3

Yeah, you apparently have the impression that you need to work on vector input and return vectors as well. But that’s not how this works. Your function should take scalar input and return scalar output.

You often see languages like matlab an python requiring vectorized functions for quadrature, for performance reasons. But that’s not optimal, adaptive quadrature/cubature functions can do better if they work on individual samples, and can also avoid array allocations.

4

It may be worth noting that Cubature.jl, as opposed to HCubature.jl, does support using functions that evaluate batches of points in vectors.

5

Hey, thanks! I changed the cop2019 function to now return return a scalar. The double integral now runs although it gives a DomainError that I’m not sure why it’s happening.

I’m rewriting all my R code in Julia in the hopes it’s a bit faster (so far I haven’t seen that happen, which is a shame) and there it works with the vectors u and v I defined. I will keep trying

6

Hi, thanks! I think I’m using Cubature.jl. Should I try with HCubature.jl instead?

7

HCubature is pure-Julia code and is a bit more flexible in some ways. e.g. it can work more types of return values, such as SVector-valued functions. Cubature, on the other hand, supports batched integrand evaluation, which can be useful for parallelization, and also has pcubature routines. (For hcubature, they use essentially the same algorithm.)

1 Like
8

You have a x^(1/alpha-2) exponentiation, i.e. a fractional power (whose value is only known at runtime). That will give a DomainError in Julia for negative real x, unless you switch to complex numbers.

For a fixed exponent, like alpha == 3, it will be much faster to hardcode the exponent. (And you should avoid globals anyway.) And x^(1/3-2) can be computed efficiently by cbrt(x)/x^2 and it will work even for negative real x.

2 Likes