[Vinecopula-commits] r143 - pkg pkg/R pkg/src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Di Sep 22 18:37:18 CEST 2015
Author: ulf
Date: 2015-09-22 18:37:17 +0200 (Tue, 22 Sep 2015)
New Revision: 143
Modified:
pkg/DESCRIPTION
pkg/R/BiCopPar2Tau.r
pkg/R/RVineGrad.r
pkg/R/RVineHessian.r
pkg/R/RVineMLE.R
pkg/R/RVinePIT.r
pkg/R/RVinePartialcorr.R
pkg/src/cdvine.c
pkg/src/deriv.c
pkg/src/evCopula.c
pkg/src/gof.c
pkg/src/hfunc.c
pkg/src/incompleteBeta.c
pkg/src/likelihood.c
pkg/src/logderiv.c
pkg/src/rvine.c
pkg/src/rvinederiv.c
pkg/src/rvinederiv2.c
pkg/src/tcopuladeriv_new.c
pkg/src/tools.c
tests/testCheck.r
tests/testMain.r
tests/testRun.r
Log:
* one new test
* DESCRIPTION file: Authors at R added
* deleted not necessary lines in code
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2015-09-21 16:16:45 UTC (rev 142)
+++ pkg/DESCRIPTION 2015-09-22 16:37:17 UTC (rev 143)
@@ -3,6 +3,20 @@
Title: Statistical Inference of Vine Copulas
Version: 1.7
Date: 2015-08-10
+Authors at R: c(
+ person("Schepsmeier", "Ulf", , "ulf.schepsmeier at tum.de", c("aut")),
+ person("Stoeber", "Jakob", , role = "aut"),
+ person("Brechmann", "Eike", "Christian", role = "aut"),
+ person("Graeler", "Benedikt", role="aut"),
+ person("Nagler", "Thomas", , "thomas.nagler at tum.de", c("aut")),
+ person("Erhardt", "Tobias", , "tobias.erhardt at tum.de", c("aut", "cre")),
+ person("Almeida", "Carlos", role="ctb"),
+ person("Min", "Aleksey", role=c("ctb", "ths")),
+ person("Czado", "Claudia", role=c("ctb", "ths")),
+ person("Hofmann", "Mathias", role="ctb"),
+ person("Killiches", "Matthias", role="ctb"),
+ person("Joe", "Harry", role="ctb")
+ )
Author: Ulf Schepsmeier, Jakob Stoeber, Eike Christian Brechmann, Benedikt Graeler, Thomas Nagler, Tobias Erhardt
Maintainer: Tobias Erhardt <tobias.erhardt at tum.de>
Depends: R (>= 2.11.0)
Modified: pkg/R/BiCopPar2Tau.r
===================================================================
--- pkg/R/BiCopPar2Tau.r 2015-09-21 16:16:45 UTC (rev 142)
+++ pkg/R/BiCopPar2Tau.r 2015-09-22 16:37:17 UTC (rev 143)
@@ -84,16 +84,14 @@
} else if (family == 9 || family == 19) {
theta <- par
delta <- par2
- # tau=1-2/(delta*(2-theta))+4/(theta^2*delta)*gamma(delta+2)*gamma((2-2*theta)/(theta)+1)/gamma(delta+3+(2-2*theta)/(theta))
+
kt <- function(t, th, de) {
((1 - (1 - t)^th)^-de - 1)/(-th * de * (1 - t)^(th - 1) * (1 - (1 - t)^th)^(-de - 1))
}
tau <- 1 + 4 * mapply(function(theta, delta) {
integrate(function(t) kt(t, th = theta, de = delta), 0, 1)$value
}, theta, delta)
- # kt <- function(t) { 1/( (1-t)^(par-1) ) } kt2 <- function(t) { 1-t } kt3 <-
- # function(t) { 1/( (1-t)^(par-1)*(1-(1-t)^par)^(-par2-1) ) } tau <-
- # 1-4/par/par2*(integrate(kt,0,1)$value-integrate(kt2,0,1)$value-integrate(kt3,0,1)$value)
+
} else if (family == 10 || family == 20) {
theta <- par
delta <- par2
@@ -110,8 +108,7 @@
} else if (family == 24 || family == 34) {
tau <- -1 - 1/par
} else if (family == 26 || family == 36) {
- # tau <- -1-4/par^2*integrate(function(x) log(x)*x*(1-x)^(2*(1+par)/-par), 0,
- # 1)$value
+
theta <- -par
param1 <- 2/theta + 1
tem <- digamma(2) - digamma(param1)
@@ -138,8 +135,7 @@
} else if (family == 29 || family == 39) {
theta <- -par
delta <- -par2
- # tau <-
- # 1-2/(delta*(2-theta))+4/(theta^2*delta)*gamma(delta+2)*gamma((2-2*theta)/(theta)+1)/gamma(delta+3+(2-2*theta)/(theta))
+
kt <- function(t, th, de) {
((1 - (1 - t)^th)^(-de) - 1)/(-th * de * (1 - t)^(th - 1) * (1 - (1 - t)^th)^(-de - 1))
}
Modified: pkg/R/RVineGrad.r
===================================================================
--- pkg/R/RVineGrad.r 2015-09-21 16:16:45 UTC (rev 142)
+++ pkg/R/RVineGrad.r 2015-09-22 16:37:17 UTC (rev 143)
@@ -39,14 +39,12 @@
stop("'RVM' has to be an RVineMatrix object.")
- # if(any(is.na(calcupdate))) { n=dim(RVM) calcupdate=array(0,dim=c(n,n,n,n)) for(i in (n-1):1){ for(k in n:(i+1)){ calcupdate[, ,k,i
- # ]=RVineMatrixUpdate(RVM,k,i) } } }
+
o <- diag(RVM$Matrix)
if (any(o != length(o):1)) {
oldRVM <- RVM
RVM <- normalizeRVineMatrix(RVM)
- # RVM = getFromNamespace('normalizeRVineMatrix','VineCopula')(RVM)
data <- data[, o[length(o):1]]
}
@@ -65,7 +63,6 @@
ll <- as.vector(V$value)
vv <- as.vector(V$direct)
vv2 <- as.vector(V$indirect)
- # calcup=as.vector(calcupdate)
w1 <- as.vector(RVM$family)
w1[is.na(w1)] <- 0
@@ -81,7 +78,6 @@
maxmat[is.na(maxmat)] <- 0
condirect[is.na(condirect)] <- 0
conindirect[is.na(conindirect)] <- 0
- # tilde_vdirect_array=array(0,dim=c(n,n,N,n,n)) tilde_vindirect_array=array(0,dim=c(n,n,N,n,n)) tilde_value_array=array(0,dim=c(n,n,N,n,n))
out <- rep(0, sum(posParams[lower.tri(posParams, diag = FALSE)]) + sum(w1 == 2))
@@ -101,11 +97,7 @@
as.double(ll),
as.double(vv),
as.double(vv2),
- #as.integer(calcup),
as.integer(as.vector(posParams)),
- #as.double(as.vector(tilde_vdirect_array)),
- #as.double(as.vector(tilde_vindirect_array)),
- #as.double(as.vector(tilde_value_array)),
PACKAGE = 'VineCopula')
@@ -122,10 +114,7 @@
gradient <- c(gradient, grad2[tt:1])
}
- # tilde_vdirect=out[[16]] tilde_vindirect=out[[17]] tilde_value=out[[18]] V$tilde_direct = array(tilde_vdirect,dim=c(n,n,N,n,n)) V$tilde_indirect =
- # array(tilde_vindirect,dim=c(n,n,N,n,n)) V$tilde_value = array(tilde_value,dim=c(n,n,N,n,n))
-
- # out2=list(gradient=gradient,V=V)
+
out2 <- list(gradient = gradient)
return(out2)
}
Modified: pkg/R/RVineHessian.r
===================================================================
--- pkg/R/RVineHessian.r 2015-09-21 16:16:45 UTC (rev 142)
+++ pkg/R/RVineHessian.r 2015-09-22 16:37:17 UTC (rev 143)
@@ -28,7 +28,6 @@
o <- diag(RVM$Matrix)
if (any(o != length(o):1)) {
oldRVM <- RVM
- # RVM = normalizeRVineMatrix(RVM)
RVM <- getFromNamespace("normalizeRVineMatrix", "VineCopula")(RVM)
data <- data[, o[length(o):1]]
}
@@ -62,10 +61,6 @@
der <- matrix(out[[13]], dd + tt, dd + tt)
subder <- matrix(out[[14]], dd + tt, dd + tt)
- # der[1:dd,1:dd]=der[dd:1,dd:1] if(tt>0) { der[(dd+1):(dd+tt),1:dd]=der[(dd+tt):(dd+1),dd:1] der[1:dd,(dd+1):(dd+tt)]=der[dd:1,(dd+tt):(dd+1)]
- # der[(dd+1):(dd+tt),(dd+1):(dd+tt)]=der[(dd+tt):(dd+1),(dd+tt):(dd+1)] } hessian[1:dd,1:dd]=hessian[dd:1,dd:1] if(tt>0) {
- # hessian[(dd+1):(dd+tt),1:dd]=hessian[(dd+tt):(dd+1),dd:1] hessian[1:dd,(dd+1):(dd+tt)]=hessian[dd:1,(dd+tt):(dd+1)]
- # hessian[(dd+1):(dd+tt),(dd+1):(dd+tt)]=hessian[(dd+tt):(dd+1),(dd+tt):(dd+1)] }
test <- apply(hessian, 2, function(x) max(abs(x)))
hessian <- hessian[test > 0, test > 0]
Modified: pkg/R/RVineMLE.R
===================================================================
--- pkg/R/RVineMLE.R 2015-09-21 16:16:45 UTC (rev 142)
+++ pkg/R/RVineMLE.R 2015-09-22 16:37:17 UTC (rev 143)
@@ -285,14 +285,10 @@
## log-likelihood function to be maximized
optim_LL <- function(parm, data, posParams, posParams2, Copula.Types, start, start2, RVM, calcupdate = NA) {
- # calcupdate=NA V=NA
nParams <- sum(posParams, na.rm = TRUE)
nParams2 <- sum(posParams2, na.rm = TRUE)
- # for(k in 1:nParams) { if(Copula.Types[k]==0) { if(k==1){parm=c(0,parm)} else if(k>length(parm)){parm=c(parm,0)} else
- # {parm=c(parm[1:(k-1)],0,parm[k:length(parm)])} } }
-
matrixParams <- start
matrixParams2 <- start2
@@ -304,7 +300,6 @@
ll <- RVineLogLik(data, RVM, par = matrixParams, par2 = matrixParams2)
- # V=ll$V
if (is.finite(ll$loglik)) {
return(ll$loglik)
@@ -322,14 +317,6 @@
nParams <- sum(posParams, na.rm = TRUE)
nParams2 <- sum(posParams2, na.rm = TRUE)
- # outparm=parm for (i in 1:length(parm)) { handle_parm=parm handle_parm[i]=handle_parm[i]+0.000001 handle_parm2=parm
- # handle_parm2[i]=handle_parm2[i]-0.000001
- # outparm[i]=(optim_LL(handle_parm,data,posParams,posParams2,Copula.Types,start,start2,RVM,calcupdate=NA)-optim_LL(handle_parm2,data,posParams,posParams2,Copula.Types,start,start2,RVM,calcupdate=NA))/(2*0.000001)
- # } print('finite differences:') print(outparm)
-
- # for(k in 1:nParams) { if(Copula.Types[k]==0) { if(k==1){parm=c(0,parm)} else if(k>length(parm)){parm=c(parm,0)} else
- # {parm=c(parm[1:(k-1)],0,parm[k:length(parm)])} } }
-
matrixParams <- start
matrixParams2 <- start2
@@ -370,7 +357,6 @@
## optimization
if (all(Copula.Types %in% c(0, 1, 2, 3:6, 13, 14, 16, 23, 24, 26, 33, 34, 36, 43, 44)) && grad == TRUE) {
- # n=dim(RVM) calcupdate=array(0,dim=c(n,n,n,n)) for(i in (n-1):1){ for(k in n:(i+1)){ calcupdate[, ,k,i ]=RVineMatrixUpdate(RVM,k,i) } }
if (hessian == TRUE || se == TRUE) {
out1 <- optim(par = startpar,
fn = optim_LL,
Modified: pkg/R/RVinePIT.r
===================================================================
--- pkg/R/RVinePIT.r 2015-09-21 16:16:45 UTC (rev 142)
+++ pkg/R/RVinePIT.r 2015-09-22 16:37:17 UTC (rev 143)
@@ -19,34 +19,7 @@
if (is(RVM)[1] != "RVineMatrix")
stop("'RVM' has to be an RVineMatrix object.")
- #if(type=="CVine") type=1
- #else if(type=="DVine") type=2
- #else if(type=="RVine") type=0
- #if(!(type %in% c(0,1,2)) ) stop("Vine type not implemented.")
-
- #if(type==1 || type==2)
- #{
- # if(type==1)
- # vine=R2CVine(RVM)
- # else if(type==2)
- # vine=R2DVine(RVM)
-
- # tmp = .C("pit",
- # as.integer(T),
- # as.integer(d),
- # as.integer(vine$family),
- # as.integer(type),
- # as.double(vine$par),
- # as.double(vine$par2),
- # as.double(data),
- # as.double(rep(0,T*d)),
- # PACKAGE='VineCopula')[[8]]
- #
- # U <- matrix(tmp,ncol=d)
- #}
- #else {
-
o <- diag(RVM$Matrix)
if (any(o != length(o):1)) {
oldRVM <- RVM
Modified: pkg/R/RVinePartialcorr.R
===================================================================
--- pkg/R/RVinePartialcorr.R 2015-09-21 16:16:45 UTC (rev 142)
+++ pkg/R/RVinePartialcorr.R 2015-09-22 16:37:17 UTC (rev 143)
@@ -199,96 +199,3 @@
}
-#######################################
-# for immeddeate testing run as well ##
-#######################################
-
-# normalizeRVineMatrix = function(RVM){
-#
-# oldOrder = diag(RVM$Matrix)
-# Matrix = reorderRVineMatrix(RVM$Matrix)
-#
-# names <- RVM$names
-# if(is.null(names))
-# names <- paste("V",1:nrow(RVM$Matrix),sep="")
-#
-# return(RVineMatrix(Matrix, RVM$family, RVM$par, RVM$par2, names = rev(names[oldOrder])))
-# }
-#
-# reorderRVineMatrix = function(Matrix){
-# oldOrder = diag(Matrix)
-#
-# O = apply(t(1:nrow(Matrix)),2,"==", Matrix)
-#
-# for(i in 1:nrow(Matrix)){
-# Matrix[O[,oldOrder[i]]] = nrow(Matrix)-i+1
-# }
-#
-# return(Matrix)
-# }
-#
-# # examples/test cases
-# ######################
-#
-# corMat <- matrix(c(1.00, 0.17, 0.15, 0.14, 0.13,
-# 0.17, 1.00, 0.30, 0.28, 0.05,
-# 0.15, 0.30, 1.00, 0.17, 0.05,
-# 0.14, 0.28, 0.17, 1.00, 0.04,
-# 0.13, 0.05, 0.05, 0.04, 1.00),5,5)
-#
-# Matrix = matrix(c(5,2,3,1,4,
-# 0,2,3,4,1,
-# 0,0,3,4,1,
-# 0,0,0,4,1,
-# 0,0,0,0,1),5,5)
-# family = matrix(1,5,5)
-#
-# par = matrix(c(0,0.2,0.9,0.5,0.8,
-# 0, 0,0.1,0.6,0.9,
-# 0, 0, 0,0.7,0.5,
-# 0, 0, 0, 0,0.8,
-# 0, 0, 0, 0, 0),5,5)
-#
-# # define RVineMatrix object
-# RVM = RVineMatrix(Matrix,family,par)
-#
-# # adjust the un-ordered RVine
-# newRVM <- RVineCor2pcor(RVM, corMat)
-# round(cor(qnorm(RVineSim(100000, newRVM)))-corMat, 2)
-#
-# # normalise the RVine
-# normRVM <- normalizeRVineMatrix(RVM)
-#
-# # adjust the normalised RVine
-# newNormRVM <- RVineCor2pcor(normRVM, corMat)
-#
-# # newRVM and newNormRVM are the same vine using only different naming:
-# newNormRVM$par - newRVM$par
-#
-# # the variable now do have a different ordering in the correlation matrix
-# newNormCor <- cor(qnorm(RVineSim(100000, newNormRVM)))
-# round(newNormCor,2)
-#
-# # permuted, they meet the initial correlation matrix up to +/- 0.01
-# round(newNormCor[c(1,4,3,2,5),c(1,4,3,2,5)]-corMat, 2)
-#
-# # re-order names of the normalised RVine generating a new RVine
-# normRVM2 <- normRVM
-# normRVM2$names <- c("V1", "V2", "V3", "V4", "V5")
-#
-# # adjust the normalised RVine
-# newNormRVM2 <- RVineCor2pcor(normRVM2, corMat)
-# # check whether the parameters are different beyond permutation (that's why
-# # permutation does not work)
-# newNormRVM2$par
-# newRVM$par
-#
-# # adjust the normalised RVine
-# newNormRVM2 <- RVineCor2pcor(normRVM2, corMat[c(1,4,3,2,5),c(1,4,3,2,5)])
-# # check whether the parameters are now identical
-# round(newNormRVM2$par - newRVM$par,2)
-#
-# # back and forth
-# RVinePcor2cor(RVineCor2pcor(RVM, corMat))-corMat
-# RVinePcor2cor(RVineCor2pcor(normRVM, corMat))-corMat
-# RVinePcor2cor(RVineCor2pcor(normRVM2, corMat))-corMat
Modified: pkg/src/cdvine.c
===================================================================
--- pkg/src/cdvine.c 2015-09-21 16:16:45 UTC (rev 142)
+++ pkg/src/cdvine.c 2015-09-22 16:37:17 UTC (rev 143)
@@ -1,611 +1,603 @@
-/*
-** cdvine.c - C code of the package CDRVine
-**
-** with contributions from Carlos Almeida, Aleksey Min,
-** Ulf Schepsmeier, Jakob Stoeber and Eike Brechmann
-**
-** A first version was based on code
-** from Daniel Berg <daniel at danielberg.no>
-** provided by personal communication.
-**
-*/
-
-
-// Include all the head files
-#include "include/vine.h" // general one
-#include "include/memoryhandling.h" // for creating two and three dimensional arrays
-#include "include/likelihood.h" // formally main functionality; log-likelihood with help functions; bivariate densities
-#include "include/cdvine.h" // Header file for this C-file
-#include "include/hfunc.h" // h-functions, i.e. conditional densities; also inverse h-functions
-
-#define UMAX 1-1e-10
-
-#define UMIN 1e-10
-
-#define XEPS 1e-4
-
-
-//////////////////////////////////////////////////////////////
-// Function to simulate from a C- or D-vine
-// Input:
-// n sample size
-// d dimension (>= 2)
-// type vine type (1=Canonical vine, 2=D-vine)
-// family copula family (see help pages which families are now included)
-// par parameter values (at least d*(d-1)/2 parameters)
-// nu second parameter for t-copula, BB-copulas and Tawn
-//
-// Output:
-// out two dimensional array of simulated data
-////////////////////////////////////////////////////////////////
-
-void pcc(int* n, int* d, int* family, int* type, double* par, double* nu, double* out)
-{
- int i, j, in=1, k, **fam;
- double *w, **v, t, **theta, **x, **ny;
-
- GetRNGstate(); //Init random number generator
- //Allocate memory:
- w = Calloc((*d+1),double);
-
- v = create_matrix(*d+1,2*(*d)-1);
- theta = create_matrix(*d,*d);
- x = create_matrix(*n+1,*d+1);
- ny = create_matrix(*d,*d);
- fam = create_intmatrix(*d,*d);
- //Initialize dependency parameters
-
- // The function arguments are one-dimensional vectors; for better understanding the transform them back to matrices (see theory)
- // This step may be updated in the future to optimize the algorithms
- k = 0;
- for(i=1;i<=*d-1;i++)
- {
- for(j=1;j<=*d-i;j++)
- {
- fam[i][j] = family[k];
- ny[i][j] = nu[k];
- theta[i][j] = par[k];
- k ++;
- }
- }
- //Simulate: (it follows the theoretical algorithm)
- if(*type==1) //Canonical vine
- {
- for(j=1;j<=*n;j++) // run over all observations (rows)
- {
- for(i=1;i<=*d;i++) w[i] = runif(0,1);
- x[j][1] = w[1];
- for(i=2;i<=*d;i++) // run over all dimensions (cols)
- {
- t = w[i];
- for(k=i-1;k>=1;k--)
- {
- //Hinv1(&fam[k][i-k],&in, &v[i][1],&v[k][k],&theta[k][i-k],&ny[k][i-k],&v[i][1]);
- Hinv1(&fam[k][i-k],&in, &t,&w[k],&theta[k][i-k],&ny[k][i-k],&t);
- }
- x[j][i] = t;
- /*if(i<*d)
- {
- for(k=1;k<i;k++) Hfunc1(&fam[k][i-k],&in, &v[i][k],&v[k][k],&theta[k][i-k],&ny[k][i-k],&v[i][k+1]);
- } */
- }
- }
- }
- else if(*type==2) //D-vine
- {
- for(j=1;j<=*n;j++)
- {
- for(i=1;i<=*d;i++) { w[i] = runif(0,1);} // printf("%10.8f \t",w[i]);} ; printf("\n");
- v[1][1] = w[1];
- v[2][1] = w[2];
- //printf("inv: %d,%d : %5.2f : %10.8f \t",1,0, theta[1][1], v[1][1]);
- Hinv1(&fam[1][1],&in,&w[2],&v[1][1],&theta[1][1],&ny[1][1],&v[2][1]);
- //printf("%10.8f \n",v[2][1]);
- Hfunc2(&fam[1][1],&in, &v[1][1],&v[2][1],&theta[1][1],&ny[1][1],&v[2][2]);
- for(i=3;i<=*d;i++)
- {
- v[i][1] = w[i];
-
- for(k=i-1;k>=2;k--)
- {
- //printf("inv: %d,%d : %5.2f : %10.8f \t",k,i-k, theta[k][i-k], v[i-1][2*k-2]);
- Hinv1(&fam[k][i-k],&in, &v[i][1],&v[i-1][2*k-2],&theta[k][i-k],&ny[k][i-k],&v[i][1]);
- //printf("%10.8f \n",v[i][1]);
- }
- Hinv1(&fam[1][i-1],&in, &v[i][1],&v[i-1][1],&theta[1][i-1],&ny[1][i-1],&v[i][1]);
- //printf("inv: %d,%d : %5.2f : %10.8f \t %10.8f \n",1,i-1, theta[1][i-1], v[i-1][1],v[i][1]);
- // Compute conditional cdf's needed in next step:
- if(i<*d)
- {
- Hfunc2(&fam[1][i-1],&in, &v[i-1][1],&v[i][1],&theta[1][i-1],&ny[1][i-1],&v[i][2]);
- Hfunc1(&fam[1][i-1],&in, &v[i][1],&v[i-1][1],&theta[1][i-1],&ny[1][i-1],&v[i][3]);
- if(i>3)
- {
- for(k=2;k<=(i-2);k++)
- {
- Hfunc2(&fam[k][i-k],&in, &v[i-1][2*k-2],&v[i][2*k-1],&theta[k][i-k],&ny[k][i-k],&v[i][2*k]);
- Hfunc1(&fam[k][i-k],&in, &v[i][2*k-1],&v[i-1][2*k-2],&theta[k][i-k],&ny[k][i-k],&v[i][2*k+1]);
- }
- }
- Hfunc2(&fam[i-1][1],&in, &v[i-1][2*i-4],&v[i][2*i-3],&theta[i-1][1],&ny[i-1][1],&v[i][2*i-2]);
- }
- }
- for(i=1;i<=*d;i++) x[j][i] = v[i][1];
- }
- }
- //Write to output vector:
- k = 0;
- for(i=1;i<=*d;i++)
- {
- for(j=1;j<=*n;j++)
- {
- out[k] = x[j][i];
- k ++;
- }
- }
- PutRNGstate(); // Function for the random number generator
- //Free memory:
- Free(w); free_matrix(v,*d+1); free_matrix(theta,*d); free_matrix(ny,*d); free_intmatrix(fam,*d); free_matrix(x,*n+1);
-}
-
-
-
-
-//////////////////////////////////////////////////////////////
-// Function to compute -log-likelihood for C- and D-vine
-// Input:
-// n sample size
-// d dimension (>=2)
-// type vine type (1=canonical vine, 2=d-vine)
-// family copula families
-// par parameter values (at least d*(d-1)/2 parameters ) The second parameter is added at the end of par
-// data data set for which to compute log-likelihood
-// Output:
-// out Log-likelihood
-// ll array with the contribution to LL (for each copula)
-// vv array for the transformation operated (Hfunc)
-/////////////////////////////////////////////////////////////
-void VineLogLikm(int* T, int* d, int* type, int* family, double* par, double* data,
- double* out, double* ll, double* vv)
-{
- int i, j, k, t, kk, **fam;
- double loglik=0.0, sumloglik=0.0, **x, **theta, **nu, ***v;
-
- //Allocate memory:
- x = create_matrix(*d+1,*T);
- // By Ulf Schepsmeier
- if(*type==1) //C-vine
- {
- v = create_3darray(*d-1,*d,*T);
- }
- else //D-vine
- {
- v = create_3darray(*d,2*(*d)-3,*T);
- }
-
- theta = create_matrix(*d,*d);
- nu = create_matrix(*d,*d);
- fam = create_intmatrix(*d+1,*d+1);
-
- //Initialize:
- k = 0;
- for(i=1;i<=*d;i++)
- {
- for (t=0;t<=*T-1;t++ )
- {
- x[i][t] = data[k]; //transform the data back into a 2-dim array
- k++;
- }
- }
- k = 0;
- for(i=1;i<=(*d-1);i++)
- {
- for(j=1;j<=(*d-i);j++)
- {
- theta[i][j] = par[k];
- fam[i][j] = family[k];
- nu[i][j] = par[*d*(*d-1)/2+k]; // the second parameter is added at the end of par (not the best solution but was practise at the beginning)
- k++;
- }
- }
-
- if(*type==1) //C-vine
- {
- // By Ulf Schepsmeier
- kk=0;
- //Compute likelihood at level 1:
- for(i=1;i<*d;i++)
- {
- LL_mod2(&fam[1][i],T,x[1],x[i+1],&theta[1][i],&nu[1][i],&loglik); // call the bivariate log-likelihood function
- //(with the correct rotation for 90, 180 and 270 degrees)
- sumloglik += loglik; // sum up
- ll[kk] = loglik; // store all bivariate log-likelihoods too
- ++kk;
- if(*d>2)
- {
- //Compute variables for next level:
- Hfunc1(&fam[1][i],T,x[i+1],x[1],&theta[1][i],&nu[1][i],v[1][i]);
- }
- }
- //Compute likelihood at next levels:
- if(*d>2)
- {
- for(k=2;k<=(*d-1);k++)
- {
- for(i=1;i<=(*d-k);i++)
- {
- LL_mod2(&fam[k][i],T,v[k-1][1],v[k-1][i+1],&theta[k][i],&nu[k][i],&loglik);
- sumloglik += loglik;
- ll[kk] = loglik;
- ++kk;
- }
- if(k<(*d-1))
- {
- for (i=1;i<=(*d-k);i++)
- {
- Hfunc1(&fam[k][i],T,v[k-1][i+1],v[k-1][1],&theta[k][i],&nu[k][i],v[k][i]);
- }
- }
- }
- }
- }
- else if(*type==2) //D-vine
- {
- kk=0;
- //Compute the likelihood at level 1:
- for(i=1;i<*d;i++)
- {
- LL_mod2(&fam[1][i],T,x[i],x[i+1],&theta[1][i],&nu[1][i],&loglik);
- sumloglik += loglik;
- ll[kk] = loglik;
- ++kk;
- }
- //Compute variables for next level:
- if(*d>2)
- {
- Hfunc2(&fam[1][1],T,x[1],x[2],&theta[1][1],&nu[1][1],v[1][1]);
- for(k=1;k<=(*d-3);k++)
- {
- Hfunc1(&fam[1][k+1],T,x[k+2],x[k+1],&theta[1][k+1],&nu[1][k+1],v[1][2*k]);
- Hfunc2(&fam[1][k+1],T,x[k+1],x[k+2],&theta[1][k+1],&nu[1][k+1],v[1][2*k+1]);
- }
- Hfunc1(&fam[1][*d-1],T,x[*d],x[*d-1],&theta[1][*d-1],&nu[1][*d-1],v[1][2*(*d)-4]);
- //Compute likelihood at next levels:
- for(k=2;k<=(*d-1);k++)
- {
- for(i=1;i<=(*d-k);i++)
- {
- LL_mod2(&fam[k][i],T,v[k-1][2*i-1],v[k-1][2*i],&theta[k][i],&nu[k][i],&loglik);
- sumloglik += loglik;
- ll[kk] = loglik; ++kk;
- }
- if(k<(*d-1))
- {
- Hfunc2(&fam[k][1],T,v[k-1][1],v[k-1][2],&theta[k][1],&nu[k][1],v[k][1]);
- if((*d)>4)
- {
- for(i=1;i<=(*d-k-2);i++)
- {
- Hfunc1(&fam[k][i+1],T,v[k-1][2*i+2],v[k-1][2*i+1],&theta[k][i+1],&nu[k][i+1],v[k][2*i]);
- Hfunc2(&fam[k][i+1],T,v[k-1][2*i+1],v[k-1][2*i+2],&theta[k][i+1],&nu[k][i+1],v[k][2*i+1]);
- }
- }
- Hfunc1(&fam[k][*d-k],T,v[k-1][2*(*d)-2*k],v[k-1][2*(*d)-2*k-1],&theta[k][*d-k],&nu[k][*d-k],v[k][2*(*d)-2*k-2]);
- }
- }
- }
- }
- //Write to output:
- *out = -sumloglik;
- kk=00;
- // By Ulf Schepsmeier
- if(*type==1) //C-Vine
- {
- if(*d>2)
- {
- for(k=1;k<(*d-1);k++)
- for(i=1;i<=(*d-k);i++)
- for(t=0;t<*T;t++)
- {
- vv[kk] = v[k][i][t]; // transformation from a 3-dim array to a vector
- ++kk;
- }
- }
- //Free memory:
- free_3darray(v,*d-1,*d);
- }
- else //D-Vine
- {
- if(*d>2)
- {
- for(k=1;k<*d;k++)
- for(i=1;i<=2*(*d-k-1);i++)
- for(t=0;t<*T;t++)
- {
- vv[kk] = v[k][i][t];
- ++kk;
- }
- }
- //Free memory:
- free_3darray(v,*d,2*(*d)-3);
- }
- //Free memory:
- free_matrix(x,*d+1); free_matrix(theta,*d); free_matrix(nu,*d); free_intmatrix(fam,*d+1);
-}
-
-//////////////////////////////////////////////////////////////
-// Function to compute an update of the log-likelihood for C- and D-vine
-// Input:
-// n sample size
-// d dimension (>=2)
-// type vine type (1=canonical vine, 2=d-vine)
-// family copula families
-// par parameter values (at least d*(d-1)/2 parameters
-// mpar index of modified parameter (related to previous computation)
-// data data set for which to compute log-likelihood
-// ll array with the stored contribution of the likelihood in a previous computation
-// vv 3d array array with the stored transformations in a previous computation
-// Output:
-// out log-likelihood (updated)
-// ll array with the contribution to LL (for each copula)
-// vv array for the transformation operated (Hfunc)
-/////////////////////////////////////////////////////////////
-void VineLogLikmP(int* T, int* d, int* type, int* family, double* par, int* mpar, double* data,
- double* out, double* ll, double* vv)
-{
- int i, j, ii, jj, k, t, kk,**fam, **ind;
- double sumloglik=0.0,loglik=0.0, **x, **theta, **nu, ***v;
- //Allocate memory:
- x = create_matrix(*d+1,*T);
-
- // By Ulf Schepsmeier
- if(*type==1) //C-vine
- {
- v = create_3darray(*d-1,*d,*T);
- }
- else //D-vine
- {
- v = create_3darray(*d,2*(*d)-3,*T);
- }
-
- theta = create_matrix(*d,*d);
- nu = create_matrix(*d,*d);
- fam = create_intmatrix(*d,*d);
- ind = create_intmatrix(*d,*d);
- //Initialize:
- k = 0;
- for(i=1;i<=*d;i++)
- {
- for (t=0;t<=*T-1;t++ )
- {
- x[i][t] = data[k];
- k++;
- }
- }
- k = 0;
- jj = *d;
- ii = *d;
- kk=00;
-
- // By Ulf Schepsmeier
- if(*type==1) //C-Vine
- {
- for(i=1;i<=(*d-1);i++)
- {
- for(j=1;j<=(*d-1);j++)
- {
- ind[i][j] = 0;
- if(j <= *d-i)
- {
- ++k;
- if (k == *mpar)
- {
- ii=i;
- jj=j;
- }
- if(ii+jj-i>0)
- {
- if(i>=ii && j==ii+jj-i)
- {
- ind[i][j]=1;
- }
- }
- else if(ii+jj-i<=0)
- {
- ind[i][j]=1;
- }
- }
- }
- }
-
- for(k=1;k<*d-1;k++)
- for(i=1;i<=*d-k;i++)
- for(t=0;t<*T;t++)
- {
- v[k][i][t] = vv[kk];
- ++kk;
- }
- }
- else //D-Vine
- {
- for(i=1;i<=(*d-1);i++)
- {
- for(j=1;j<=(*d-1);j++)
- {
- ind[i][j] = 0;
- if(j <= *d-i)
- {
- ++k;
- if (k == *mpar)
- {
- ii=i;
- jj=j;
- }
- if(i >= ii && j >= ii+jj-i && j <= *d-i && j <= jj)
- {
- ind[i][j]=1;
- }
- // //printf("%d ", ind[i][j]);
- }
- }
- // //printf("\n");
- }
-
- for(k=1;k<*d;k++)
- for(i=1;i<=2*(*d-k-1);i++)
- for(t=0;t<*T;t++)
- {
- v[k][i][t] = vv[kk];
- ++kk;
- }
- }
- k=0;
- for(i=1;i<=(*d-1);i++)
- {
- for(j=1;j<=(*d-i);j++)
- {
- theta[i][j] = par[k];
- fam[i][j] = family[k];
- nu[i][j] = par[*d*(*d-1)/2+k];
- k ++;
- }
- }
-
-
-
- if(*type==1) //C-vine
- {
- // By Ulf Schepsmeier
- kk=0;
- //Compute likelihood at level 1:
- for(i=1;i<*d;i++)
- {
- if(ind[1][i]==1)
- {
- LL_mod2(&fam[1][i],T,x[1],x[i+1],&theta[1][i],&nu[1][i],&loglik);
- ll[kk] = loglik;
- if(*d>2)
- {
- //Compute variables for next level:
- Hfunc1(&fam[1][i],T,x[i+1],x[1],&theta[1][i],&nu[1][i],v[1][i]);
- }
- }
- sumloglik += ll[kk];
- ++kk;
- }
- //Compute likelihood at next levels:
- if(*d>2)
- {
- for(k=2;k<=(*d-1);k++)
- {
- for(i=1;i<=(*d-k);i++)
- {
- if(ind[k][i]==1)
- {
- LL_mod2(&fam[k][i],T,v[k-1][1],v[k-1][i+1],&theta[k][i],&nu[k][i],&loglik);
- ll[kk] = loglik;
- if(k<(*d-1))
- {
- Hfunc1(&fam[k][i],T,v[k-1][i+1],v[k-1][1],&theta[k][i],&nu[k][i],v[k][i]);
- }
- }
- sumloglik += ll[kk];
- ++kk;
- }
- }
- }
- }
- else if(*type==2) //D-vine
- {
- kk=0;
- //Compute the likelihood at level 1:
- for(i=1;i<*d;i++)
- {
- if(ind[1][i]==1)
- {
- LL_mod2(&fam[1][i],T,x[i],x[i+1],&theta[1][i],&nu[1][i],&loglik);
- ll[kk] = loglik;
- }
- sumloglik += ll[kk];
- ++kk;
- }
- //Compute variables for next level:
- if(*d>2)
- {
- if(ind[1][1]==1) Hfunc2(&fam[1][1],T,x[1],x[2],&theta[1][1],&nu[1][1],v[1][1]);
- for(k=1;k<=(*d-3);k++)
- if(ind[1][k+1]==1)
- {
- Hfunc1(&fam[1][k+1],T,x[k+2],x[k+1],&theta[1][k+1],&nu[1][k+1],v[1][2*k]);
- Hfunc2(&fam[1][k+1],T,x[k+1],x[k+2],&theta[1][k+1],&nu[1][k+1],v[1][2*k+1]);
- }
- if(ind[1][*d-1]) Hfunc1(&fam[1][*d-1],T,x[*d],x[*d-1],&theta[1][*d-1],&nu[1][*d-1],v[1][2*(*d)-4]);
- //Compute likelihood at next levels:
- for(k=2;k<=(*d-1);k++)
- {
- for(i=1;i<=(*d-k);i++)
- {
- if(ind[k][i] == 1)
- {
- LL_mod2(&fam[k][i],T,v[k-1][2*i-1],v[k-1][2*i],&theta[k][i],&nu[k][i],&loglik);
- ll[kk] = loglik;
- }
- sumloglik += ll[kk];
- ++kk;
- }
- if(k<(*d-1))
- {
- if(ind[k][1] == 1) Hfunc2(&fam[k][1],T,v[k-1][1],v[k-1][2],&theta[k][1],&nu[k][1],v[k][1]);
- if((*d)>4)
- {
- for(i=1;i<=(*d-k-2);i++)
- {
- if(ind[k][i+1]==1)
- {
- Hfunc1(&fam[k][i+1],T,v[k-1][2*i+2],v[k-1][2*i+1],&theta[k][i+1],&nu[k][i+1],v[k][2*i]);
- Hfunc2(&fam[k][i+1],T,v[k-1][2*i+1],v[k-1][2*i+2],&theta[k][i+1],&nu[k][i+1],v[k][2*i+1]);
- }
- }
- }
- if(ind[k][*d-k]==1) Hfunc1(&fam[k][*d-k],T,v[k-1][2*(*d)-2*k],v[k-1][2*(*d)-2*k-1],&theta[k][*d-k],&nu[k][*d-k],v[k][2*(*d)-2*k-2]);
- }
- }
- }
- }
- //Write to output:
- *out = -sumloglik;
- kk=00;
- // By Ulf Schepsmeier
- if(*type==1) //C-Vine
- {
- if(*d>2)
- {
- for(k=1;k<*d-1;k++)
- for(i=1;i<=*d-k;i++)
- for(t=0;t<*T;t++)
- {
- vv[kk] = v[k][i][t];
- ++kk;
- }
- }
- //Free memory:
- free_3darray(v,*d-1,*d);
- }
- else //D-Vine
- {
- if(*d>2)
- {
- for(k=1;k<*d;k++)
- for(i=1;i<=2*(*d-k-1);i++)
- for(t=0;t<*T;t++)
- {
- vv[kk] = v[k][i][t];
- ++kk;
- }
- }
- //Free memory:
- free_3darray(v,*d,2*(*d)-3);
- }
- //Free memory:
- free_matrix(x,*d+1); free_matrix(theta,*d); free_matrix(nu,*d); free_intmatrix(fam,*d);free_intmatrix(ind,*d);
-}
+/*
+** cdvine.c - C code of the package CDRVine
+**
+** with contributions from Carlos Almeida, Aleksey Min,
+** Ulf Schepsmeier, Jakob Stoeber and Eike Brechmann
+**
+** A first version was based on code
+** from Daniel Berg <daniel at danielberg.no>
+** provided by personal communication.
+**
+*/
+
+
+// Include all the head files
+#include "include/vine.h" // general one
+#include "include/memoryhandling.h" // for creating two and three dimensional arrays
+#include "include/likelihood.h" // formally main functionality; log-likelihood with help functions; bivariate densities
+#include "include/cdvine.h" // Header file for this C-file
+#include "include/hfunc.h" // h-functions, i.e. conditional densities; also inverse h-functions
+
+#define UMAX 1-1e-10
+
+#define UMIN 1e-10
+
+#define XEPS 1e-4
+
+
+//////////////////////////////////////////////////////////////
+// Function to simulate from a C- or D-vine
+// Input:
+// n sample size
+// d dimension (>= 2)
+// type vine type (1=Canonical vine, 2=D-vine)
+// family copula family (see help pages which families are now included)
+// par parameter values (at least d*(d-1)/2 parameters)
+// nu second parameter for t-copula, BB-copulas and Tawn
+//
+// Output:
+// out two dimensional array of simulated data
+////////////////////////////////////////////////////////////////
+
+void pcc(int* n, int* d, int* family, int* type, double* par, double* nu, double* out)
+{
+ int i, j, in=1, k, **fam;
+ double *w, **v, t, **theta, **x, **ny;
+
+ GetRNGstate(); //Init random number generator
+ //Allocate memory:
+ w = Calloc((*d+1),double);
+
+ v = create_matrix(*d+1,2*(*d)-1);
+ theta = create_matrix(*d,*d);
+ x = create_matrix(*n+1,*d+1);
+ ny = create_matrix(*d,*d);
+ fam = create_intmatrix(*d,*d);
+ //Initialize dependency parameters
+
+ // The function arguments are one-dimensional vectors; for better understanding the transform them back to matrices (see theory)
+ // This step may be updated in the future to optimize the algorithms
+ k = 0;
+ for(i=1;i<=*d-1;i++)
+ {
+ for(j=1;j<=*d-i;j++)
+ {
+ fam[i][j] = family[k];
+ ny[i][j] = nu[k];
+ theta[i][j] = par[k];
+ k ++;
+ }
+ }
+ //Simulate: (it follows the theoretical algorithm)
+ if(*type==1) //Canonical vine
+ {
+ for(j=1;j<=*n;j++) // run over all observations (rows)
+ {
+ for(i=1;i<=*d;i++) w[i] = runif(0,1);
+ x[j][1] = w[1];
+ for(i=2;i<=*d;i++) // run over all dimensions (cols)
+ {
+ t = w[i];
+ for(k=i-1;k>=1;k--)
+ {
+ Hinv1(&fam[k][i-k],&in, &t,&w[k],&theta[k][i-k],&ny[k][i-k],&t);
+ }
+ x[j][i] = t;
+ }
+ }
+ }
+ else if(*type==2) //D-vine
+ {
+ for(j=1;j<=*n;j++)
+ {
+ for(i=1;i<=*d;i++) { w[i] = runif(0,1);}
+ v[1][1] = w[1];
+ v[2][1] = w[2];
+
+ Hinv1(&fam[1][1],&in,&w[2],&v[1][1],&theta[1][1],&ny[1][1],&v[2][1]);
+ Hfunc2(&fam[1][1],&in, &v[1][1],&v[2][1],&theta[1][1],&ny[1][1],&v[2][2]);
+ for(i=3;i<=*d;i++)
+ {
+ v[i][1] = w[i];
+
+ for(k=i-1;k>=2;k--)
+ {
+ Hinv1(&fam[k][i-k],&in, &v[i][1],&v[i-1][2*k-2],&theta[k][i-k],&ny[k][i-k],&v[i][1]);
+ }
+ Hinv1(&fam[1][i-1],&in, &v[i][1],&v[i-1][1],&theta[1][i-1],&ny[1][i-1],&v[i][1]);
+
+ // Compute conditional cdf's needed in next step:
+ if(i<*d)
+ {
+ Hfunc2(&fam[1][i-1],&in, &v[i-1][1],&v[i][1],&theta[1][i-1],&ny[1][i-1],&v[i][2]);
+ Hfunc1(&fam[1][i-1],&in, &v[i][1],&v[i-1][1],&theta[1][i-1],&ny[1][i-1],&v[i][3]);
+ if(i>3)
+ {
+ for(k=2;k<=(i-2);k++)
+ {
+ Hfunc2(&fam[k][i-k],&in, &v[i-1][2*k-2],&v[i][2*k-1],&theta[k][i-k],&ny[k][i-k],&v[i][2*k]);
+ Hfunc1(&fam[k][i-k],&in, &v[i][2*k-1],&v[i-1][2*k-2],&theta[k][i-k],&ny[k][i-k],&v[i][2*k+1]);
+ }
+ }
+ Hfunc2(&fam[i-1][1],&in, &v[i-1][2*i-4],&v[i][2*i-3],&theta[i-1][1],&ny[i-1][1],&v[i][2*i-2]);
+ }
+ }
+ for(i=1;i<=*d;i++) x[j][i] = v[i][1];
+ }
+ }
+ //Write to output vector:
+ k = 0;
+ for(i=1;i<=*d;i++)
+ {
+ for(j=1;j<=*n;j++)
+ {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vinecopula -r 143
Mehr Informationen über die Mailingliste Vinecopula-commits