flr / FLa4a

The repository for the JRC initiative on stock assessment
https://fishreg.jrc.ec.europa.eu/web/a4a
12 stars 6 forks source link

mvrtriangle issue #83

Closed jaado closed 7 years ago

jaado commented 8 years ago

mvrtriangle should be a special case of mvrcop, where the copula is one of type "ellipCopula" and family "t". However, the following two copula distributions (vbObj1 and vbObj2) do not match:

Set up the object:

mm <- matrix(NA, ncol=3, nrow=3) diag(mm) <- c(50, 0.001,0.001) mm[upper.tri(mm)] <- mm[lower.tri(mm)] <- c(0.1,0.01,0.00004) vbObj <- a4aGr(grMod=~linf(1-exp(-k(t-t0))), grInvMod=~t0-1/k*log(1-len/linf), params=FLPar(linf=58.5, k=0.086, t0=0.001, units=c("cm","yr^-1","yr")), vcov=mm, distr="norm") pars <- list(list(a=50, b=100, c=58.5), list(a=0.06, b=0.2, c=0.086), list(a=0, b=0.005, c=0.001))

Apply mvrtriangle...

set.seed(1) vbObj1 <- mvrtriangle(10000, vbObj, paramMargins=pars) splom(data.frame(t(params(vbObj1)@.Data)), pch=".")

...and compare with mvrcop

set.seed(1) vbObj2 <- mvrcop(10000, vbObj, copula="ellipCopula", family="t", param=0, margins="triangle", paramMargins=pars) splom(data.frame(t(params(vbObj2)@.Data)), pch=".")

ejardim commented 7 years ago

Seems that triangle is scaling the vcov and using a uniform while rcop is not.

ejardim commented 7 years ago

To generate the same copulas one has to pass the arguments about the copula covariance.

set.seed(1) vbObj1 <- mvrtriangle(10000, vbObj, paramMargins=pars, dispstr="ex", param=0) set.seed(1) vbObj2 <- mvrcop(10000, vbObj, copula="ellipCopula", family="t", param=0, margins="triangle", paramMargins=pars) all.equal(vbObj2, vbObj1)