I did a quick check on the updated `cdf()`

and was puzzled by the finding. The result of `cdf()`

on the Gaussian copula is inconsistent with results from three alternative ways of computing the CDF. Meanwhile, the three alternatives are consistent among themselves.

```
using Copulas, BivariateCopulas, Distributions, StatsBase, MvNormalCDF, Cubature
rho = 0.5
ub1, ub2 = -0.1, 0.1
x = Normal(0, 1)
y = Normal(0, 2)
#### Copulas.jl
C1 = Copulas.GaussianCopula([1 rho; rho 1])
D1 = SklarDist(C1, (x,y))
cdf(D1, [ub1, ub2]) # 0.24018536633697982
#### alternative 1: Cubature on D1
func(x) = pdf(D1, [x[1], x[2]])
Cubature.pcubature(func, [-15, -15], [ub1, ub2], reltol=sqrt(eps())) # (0.32190029773358525, 3.3549441003088987e-11)
#### alternative 2: MvNormalCDF on D1
Ī¼ = mean.(D1.m)
Ļ = collect(std.(D1.m))
Ī£ = cor2cov(D1.C.Ī£, Ļ)
MvNormalCDF.mvnormcdf(Ī¼, Ī£, [-Inf, -Inf], [ub1, ub2], m = 10000) # (0.3219097660568473, 1.5325442373789225e-5)
#### alternative 3: BivariateCopulas.jl
C2 = BivariateCopulas.Gaussian(rho)
D2 = C2(x, y);
cdf(D2, ub1, ub2) # 0.3219002977336174
```

As you can see, the first result of `0.24018536633697982`

is quite different from others.

I looked at the source code:

```
function Distributions.cdf(C::CT,u) where {CT<:Copula}
f(x) = pdf(C,x)
z = zeros(eltype(u),length(C))
return Cubature.pcubature(f,z,u,reltol=sqrt(eps()))[1]
end
```

The line on `z`

seems to indicate that the lower bound of the integration is a vector of 0. However, I think it should be the lower bounds of the marginals of `X`

and `Y`

. In the above example where both of the marginals are normal, the domain should be in the range of (-Inf, u).

How to fix the program? I am scratching my head, as I donāt know how to pull out the lower bounds of user-specified marginals.