I have a non-standard multivariate distribution function/copula (here, a copula in the MWE). I am trying to compute the density function by taking the cross-partial derivatives of the distribution/copula with respect to each variable. But I am getting some values that are not between 0 and 1. Here is a minimum working example using a bivariate Gumbel copula
using FiniteDiff, Plots, LinearAlgebra
Cop(x, θ) = exp(-((-log(x[1]))^(θ) + (-log(x[2]))^(θ))^(1/θ)) # bivariate gumbel copula as example
Cop_fixed(x) = Cop(x,3)
cross_partial = []
for i in 0.01:0.01:0.999, j in 0.01:0.01:0.999 # using small grid
append!(cross_partial, FiniteDiff.finite_difference_hessian(Cop_fixed,[i,j])[1,2])
end
plot(collect(0.01:0.01:0.999), collect(0.01:0.01:0.999), cross_partial, st=:surface)
xflip!(true)
However your finite diffrences are not that good :
using FiniteDiff, Plots, LinearAlgebra, Copulas, Distributions
Cop(x, θ) = exp(-((-log(x[1]))^(θ) + (-log(x[2]))^(θ))^(1/θ)) # bivariate gumbel copula as example
Cop_fixed(x) = Cop(x,3)
cross_partial = []
G = GumbelCopula(2,3)
truth = []
for i in 0.01:0.01:0.999, j in 0.01:0.01:0.999 # using small grid
append!(cross_partial, FiniteDiff.finite_difference_hessian(Cop_fixed,[i,j])[1,2])
append!(truth, pdf(G,[i,j]))
end
relative_error = (cross_partial .- truth)./truth
plot(collect(0.01:0.01:0.999), collect(0.01:0.01:0.999), relative_error, st=:surface)
xflip!(true)
FiniteDiff actually seemed to perform better - I figured that the forward mode in ForwardDiff would make it perform especially badly around the boundary… Is there a better autodiff mode to use for this kind of problem?
I think it really depends on your problem, since this gumbel is just a MWE. But my default option is to go with autodiff and not finitediff, which is why I was suprised.
To clarify, algorithmic differentiation (eg with ForwardDiff or Enzyme) makes no floating point cancellation/truncation errors. And it is usually as fast as numerical differentiation (eg with FiniteDiff) in forward mode, possibly much faster in reverse mode.
Here we definitely want forward mode cause the input dimension is small
No, not true. It can even have unbounded error. But same with finite differencing. It does tend to have a lot less error though.
It should usually be faster as the derivative code is generally simpler that the primal, and being able to SIMD the two together can generally make the derivative pushforward less than the 2x it would otherwise take.
Let me correct: it makes no floating point errors inherent to the differentiating method. Whereas the whole point of finite differences is to compute f(x+h) - f(x) / h and thus strike a balance between Taylor residual (h large) and floating point errors (h small). Autodiff has no h, that’s what I meant
So, in general, is there a preferred autodifferentiation method for best results around the boundary of the support of a function like the edges of the cdf here?
In the particular Gumbel case, this is exactly what Copulas.jl is doing right now. The derivative of Gumbel’s generator is not that hard to compute by hand, and adding it there would improve performance, but since it is not there the default is ForwardDiff.
