# Numerical computation density function from distribution with autodifferentiation
**URL:** https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302
**Category:** Numerics
**Created:** [March 7, 2024, 2:42pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302 "2024-03-07T14:42:11Z")
**Posts on this page:** 17
**Page:** 1
### Author:  [@Kieran\_Marray](https://discourse.julialang.org/u/Kieran_Marray)
#### Post date: [March 7, 2024, 2:42pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/1 "2024-03-07T14:42:11Z")
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
```julia
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)
```
A lot of the values look correct, but I get some values greater than 1… Anyone have any idea why?
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [March 7, 2024, 2:51pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/2 "2024-03-07T14:51:41Z")
Why do you expect these values \partial f / \partial x\_1 \partial x\_2 to always be positive?
---
### Author:  [@Kieran\_Marray](https://discourse.julialang.org/u/Kieran_Marray)
#### Post date: [March 7, 2024, 3:23pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/3 "2024-03-07T15:23:46Z")
Whoops - error in the MWE. Corrected here. For the given copula, the density should equal the cross-partial derivative
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [March 7, 2024, 3:25pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/4 "2024-03-07T15:25:37Z")
A density function does not have to be between 0 and 1, so if I understand your example correctly it doesn’t seem to be a problem
---
### Author:  [@Kieran\_Marray](https://discourse.julialang.org/u/Kieran_Marray)
#### Post date: [March 7, 2024, 3:32pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/5 "2024-03-07T15:32:20Z")
Ahh yes I clearly need to take a break. Area also integrates to approx 1. Thanks for correcting my stupid error…
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [March 7, 2024, 3:40pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/6 "2024-03-07T15:40:36Z")
No worries! Densities are tricky beasts
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [March 7, 2024, 3:55pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/7 "2024-03-07T15:55:31Z")
However your finite diffrences are not that good :
```julia
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)
```
with:
```julia
julia> maximum(abs.(relative_error))
0.0001512309848833021
julia> maximum(abs.(cross_partial .- truth))
0.00961047447613339
julia>
```
Why FiniteDiff and not ForwardDiff ?
---
### Author:  [@Kieran\_Marray](https://discourse.julialang.org/u/Kieran_Marray)
#### Post date: [March 7, 2024, 4:06pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/8 "2024-03-07T16:06:12Z")
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?
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [March 7, 2024, 4:10pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/9 "2024-03-07T16:10:33Z")
> [@Kieran\_Marray](#):
>
> FiniteDiff actually seemed to perform better - I figured that the forward mode in ForwardDiff would make it perform especially badly around the boundary…
“Performs better” in term of runtime or error ?
> [@Kieran\_Marray](#):
>
> 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.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [March 7, 2024, 5:18pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/10 "2024-03-07T17:18:41Z")
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
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [March 7, 2024, 5:22pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/11 "2024-03-07T17:22:25Z")
And in general for small vectors and matrices like these, you will get an insane speed boost from using StaticArrays.jl instead of normal arrays
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 7, 2024, 6:59pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/12 "2024-03-07T18:59:36Z")
> [@gdalle](#):
>
> To clarify, algorithmic differentiation (eg with ForwardDiff or Enzyme) _makes no floating point cancellation errors_.
No, not true. It can even have unbounded error. But same with finite differencing. It does tend to have a lot less error though.
> [@gdalle](#):
>
> And it is usually as fast as numerical differentiation (eg with FiniteDiff) in forward mode
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.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [March 7, 2024, 7:26pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/13 "2024-03-07T19:26:31Z")
> [@ChrisRackauckas](#):
>
> No, not true. It can even have unbounded error. But same with finite differencing. It does tend to have a lot less error though.
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
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [March 7, 2024, 7:29pm UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/14 "2024-03-07T19:29:35Z")
My formulation was influenced by the book I’m currently reading ^^ (Griewank & Walther 2008)
---
### Author:  [@Kieran\_Marray](https://discourse.julialang.org/u/Kieran_Marray)
#### Post date: [March 8, 2024, 9:43am UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/15 "2024-03-08T09:43:31Z")
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?
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [March 8, 2024, 10:36am UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/16 "2024-03-08T10:36:20Z")
I’m not sure about this specific case, but for this input dimension ForwardDiff.jl is your best choice in general
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [March 8, 2024, 10:44am UTC](https://discourse.julialang.org/t/numerical-computation-density-function-from-distribution-with-autodifferentiation/111302/17 "2024-03-08T10:44:13Z")
I second this, ForwardDiff is great.
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](https://github.com/lrnv/Copulas.jl/blob/a91a3e64aefc5403aad337c56dad776dafa33c71/src/Generator/UnivariateGenerator/GumbelGenerator.jl#L42) would improve performance, but since it is not there the default is ForwardDiff.