[Vinecopula-commits] r46 - in pkg: R man src src/include
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Do Jan 9 13:42:24 CET 2014
Author: mhofert
Date: 2014-01-09 13:42:23 +0100 (Thu, 09 Jan 2014)
New Revision: 46
Modified:
pkg/R/RVineSim.R
pkg/man/RVineSim.Rd
pkg/src/include/rvine.h
pkg/src/rvine.c
Log:
extended RVineSim() to allow for matrix of U to be passed
Modified: pkg/R/RVineSim.R
===================================================================
--- pkg/R/RVineSim.R 2013-12-20 14:09:49 UTC (rev 45)
+++ pkg/R/RVineSim.R 2014-01-09 12:42:23 UTC (rev 46)
@@ -1,12 +1,19 @@
-RVineSim<-function(N,RVM)
-{
- if(!is(RVM, "RVineMatrix")) stop("'RVM' has to be an RVineMatrix object.")
-
+RVineSim <- function(N, RVM, U=NULL)
+{
+ stopifnot(N >= 1)
+ if(!is(RVM, "RVineMatrix")) stop("'RVM' has to be an RVineMatrix object.")
+
n = dim(RVM)
-
+
o = diag(RVM$Matrix)
RVM = normalizeRVineMatrix(RVM)
+ takeU <- !is.null(U)
+ if(takeU) {
+ if(!is.matrix(U)) U <- rbind(U, deparse.level=0L)
+ if((d <- ncol(U)) < 2) stop("U should be at least bivariate") # should be an (N, n) matrix
+ }
+
matri=as.vector(RVM$Matrix)
w1=as.vector(RVM$family)
th=as.vector(RVM$par)
@@ -32,24 +39,26 @@
as.double(th),
as.double(th2),
as.double(tmp),
+ as.double(U),
+ as.integer(takeU),
PACKAGE = 'VineCopula')[[9]]
out=matrix(tmp,ncol=n)
if(!is.null(RVM$names)){
- colnames(out) = RVM$names
+ colnames(out) = RVM$names
}
out = out[,sort(o[length(o):1],index.return=TRUE)$ix]
return(out)
-}
+}
transform<-function(M)
{
n=dim(M)[1]
-
+
M.new=matrix(rep(0,n*n),n,n)
for(i in 1:n)
{
Modified: pkg/man/RVineSim.Rd
===================================================================
--- pkg/man/RVineSim.Rd 2013-12-20 14:09:49 UTC (rev 45)
+++ pkg/man/RVineSim.Rd 2014-01-09 12:42:23 UTC (rev 46)
@@ -1,5 +1,5 @@
-\name{RVineSim}
-\alias{RVineSim}
+\name{RVineSim}
+\alias{RVineSim}
\title{Simulation from an R-vine copula model}
@@ -8,12 +8,15 @@
}
\usage{
-RVineSim(N, RVM)
+RVineSim(N, RVM, U=NULL)
}
\arguments{
\item{N}{Number of d-dimensional observations to simulate.}
- \item{RVM}{An \code{\link{RVineMatrix}} object containing the information of the R-vine copula model.}
+ \item{RVM}{An \code{\link{RVineMatrix}} object containing the
+ information of the R-vine copula model.}
+ \item{U}{If not \code{\link{NULL}}, an (N,d)-matrix of U[0,1] random variates
+ to be transformed to the copula sample.}
}
\value{
Modified: pkg/src/include/rvine.h
===================================================================
--- pkg/src/include/rvine.h 2013-12-20 14:09:49 UTC (rev 45)
+++ pkg/src/include/rvine.h 2014-01-09 12:42:23 UTC (rev 46)
@@ -5,7 +5,7 @@
// PCC for R-vine
////////////////////////////////////////////
-void SimulateRVine(int* T, int* d, int* family, int* maxmat, int* matrix, int* conindirect, double* par, double* par2, double* out);
+void SimulateRVine(int* T, int* d, int* family, int* maxmat, int* matrix, int* conindirect, double* par, double* par2, double* out, double* U, int* takeU);
//////////////////////////////////////////////////////////////
@@ -20,7 +20,7 @@
// par2 second set of parameter values (f.e. for student copulas)
// data data set for which to compute log-likelihood
// matrix an RVineMatrix in vector form
-// condirect, conindirect Matrizes which tell us where we find the right values
+// condirect, conindirect Matrizes which tell us where we find the right values
// seperate Control Parameter, do we want to seperate the likelihoods for each data point?
// calcupdate matrix which tells us for which parameters we need to redo the calculations, not newly computed values are taken from ll, vv, vv2
// seperate Control Parameter, do we want to seperate the likelihoods for each data point?
@@ -28,10 +28,10 @@
// Output:
// out Loglikelihood
// ll array with the contribution to LL (for each copula)
-// vv,vv2 array for the transformation operated (Hfunc)
+// vv,vv2 array for the transformation operated (Hfunc)
/////////////////////////////////////////////////////////////
-void VineLogLikRvine(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+void VineLogLikRvine(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
double* out, double* ll, double* vv, double* vv2, int* calcupdate, int* seperate);
@@ -47,17 +47,17 @@
// par2 second set of parameter values (f.e. for student copulas)
// data data set for which to compute log-likelihood
// matrix an RVineMatrix in vector form
-// condirect, conindirect Matrizes which tell us where we find the right values
+// condirect, conindirect Matrizes which tell us where we find the right values
// seperate Control Parameter, do we want to seperate the likelihoods for each data point?
// calcupdate matrix which tells us for which parameters we need to redo the calculations, not newly computed values are taken from ll, vv, vv2
//
// Output:
// out Loglikelihood
// ll array with the contribution to LL (for each copula)
-// vv,vv2 array for the transformation operated (Hfunc)
+// vv,vv2 array for the transformation operated (Hfunc)
/////////////////////////////////////////////////////////////
-void VineLikRvine(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+void VineLikRvine(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
double* out, double* ll, double* vv, double* vv2, int* calcupdate, int* seperate);
#endif
Modified: pkg/src/rvine.c
===================================================================
--- pkg/src/rvine.c 2013-12-20 14:09:49 UTC (rev 45)
+++ pkg/src/rvine.c 2014-01-09 12:42:23 UTC (rev 46)
@@ -1,12 +1,12 @@
/*
-** rvine.c - C code of the package CDRVine
-**
-** with contributions from Carlos Almeida, Aleksey Min,
+** rvine.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.
+** provided by personal communication.
**
*/
@@ -31,38 +31,38 @@
// par2 second set of parameter values (f.e. for student copulas)
// data data set for which to compute log-likelihood
// matrix an RVineMatrix in vector form
-// condirect, conindirect Matrizes which tell us where we find the right values
+// condirect, conindirect Matrizes which tell us where we find the right values
// seperate Control Parameter, do we want to seperate the likelihoods for each data point?
// calcupdate matrix which tells us for which parameters we need to redo the calculations, not newly computed values are taken from ll, vv, vv2
// Output:
// out Loglikelihood
// ll array with the contribution to LL (for each copula)
-// vv,vv2 array for the transformation operated (Hfunc)
+// vv,vv2 array for the transformation operated (Hfunc)
/////////////////////////////////////////////////////////////
-void VineLogLikRvine(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+void VineLogLikRvine(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
double* out, double* ll, double* vv, double* vv2, int* calcupdate, int* seperate)
{
int i, j, k, t, m, **fam;
int kk=1;
double loglik=1.0, sumloglik=0.0, **x, **theta, **nu, ***vdirect, ***vindirect;
-
+
double *sumsplitlog;
sumsplitlog=(double*) Calloc(*T,double);
-
- for (t=0;t<*T;t++ )
+
+ for (t=0;t<*T;t++ )
{
sumsplitlog[t] = 0;
}
-
-
-
+
+
+
double ***value2;
value2=create_3darray(*d,*d,*T);
double **value;
value=create_matrix(*d,*d);
-
+
//Allocate memory
x = create_matrix(*d,*T);
vdirect = create_3darray(*d,*d,*T);
@@ -75,24 +75,24 @@
//cindirect=create_intmatrix(*d,*d);
//mat=create_intmatrix(*d,*d);
//calc=create_intmatrix(*d,*d);
-
-
+
+
/* m=*d;
Rprintf("%d\n\n",m);
m=*T;
Rprintf("%d\n\n",m); */
-
+
//Initialize
k=0;
for(i=0;i<(*d);i++)
{
- for (t=0;t<*T;t++ )
+ for (t=0;t<*T;t++ )
{
x[i][t] = data[k];
k++;
}
}
-
+
k=0;
for(i=0;i<(*d);i++)
{
@@ -111,67 +111,67 @@
}
//calc[i][j]=calcupdate[(i+1)+(*d)*j-1] ;
}
- }
-
+ }
+
for(i=0;i<(*d);i++)
{
for(j=0;j<(*d);j++)
- {
- for(t=0;t<*T;t++ )
+ {
+ for(t=0;t<*T;t++ )
{
vdirect[i][j][t]=vv[(i+1)+(*d)*j+(*d)*(*d)*t-1];
vindirect[i][j][t]=vv2[(i+1)+(*d)*j+(*d)*(*d)*t-1];
if(*seperate==1)
- {
+ {
value2[i][j][t]=ll[(i+1)+(*d)*j+(*d)*(*d)*t-1];
}
}
}
}
-
-
-
+
+
+
for(i=0;i<(*d);i++)
{
- for(t=0;t<*T;t++ )
+ for(t=0;t<*T;t++ )
{
vdirect[*d-1][i][t]=x[*d-1-i][t];
}
}
-
-
+
+
for(i=*d-2; i>-1; i--)
{
for(k=*d-1;k>i;k--)
- {
+ {
//if(calc[k][i]==1)
if(calcupdate[(k+1)+(*d)*i-1]==1)
{
//m=mmat[k][i];
m=maxmat[(k+1)+(*d)*i-1];
-
+
//if(m == mat[k][i])
if(m == matrix[(k+1)+(*d)*i-1])
- {
+ {
if(*seperate==1)
{
kk = 1;
- for(t=0;t<*T;t++ )
+ for(t=0;t<*T;t++ )
{
- LL_mod(&fam[k][i],&kk,&vdirect[k][i][t],&vdirect[k][(*d-m)][t],&theta[k][i],&nu[k][i],&loglik);
+ LL_mod(&fam[k][i],&kk,&vdirect[k][i][t],&vdirect[k][(*d-m)][t],&theta[k][i],&nu[k][i],&loglik);
/* sumloglik += loglik; */
value2[k][i][t]=loglik;
}
}
else
{
-
- LL_mod(&fam[k][i],T,vdirect[k][i],vdirect[k][(*d-m)],&theta[k][i],&nu[k][i],&loglik);
+
+ LL_mod(&fam[k][i],T,vdirect[k][i],vdirect[k][(*d-m)],&theta[k][i],&nu[k][i],&loglik);
/* sumloglik += loglik; */
value[k][i]=loglik;
}
- //if(cdirect[k-1][i]==1)
+ //if(cdirect[k-1][i]==1)
if(condirect[k+(*d)*i-1]==1)
{
Hfunc1(&fam[k][i],T,vdirect[k][i],vdirect[k][*d-m],&theta[k][i],&nu[k][i],vdirect[k-1][i]);
@@ -179,29 +179,29 @@
//if(cindirect[k-1][i]==1)
if(conindirect[k+(*d)*i-1]==1)
{
- Hfunc2(&fam[k][i],T,vdirect[k][(*d-m)],vdirect[k][i],&theta[k][i],&nu[k][i],vindirect[k-1][i]);
+ Hfunc2(&fam[k][i],T,vdirect[k][(*d-m)],vdirect[k][i],&theta[k][i],&nu[k][i],vindirect[k-1][i]);
}
-
+
}
else
{
if(*seperate==1)
{
kk = 1;
- for(t=0;t<*T;t++ )
+ for(t=0;t<*T;t++ )
{
- LL_mod(&fam[k][i],&kk,&vdirect[k][i][t],&vindirect[k][(*d-m)][t],&theta[k][i],&nu[k][i],&loglik);
+ LL_mod(&fam[k][i],&kk,&vdirect[k][i][t],&vindirect[k][(*d-m)][t],&theta[k][i],&nu[k][i],&loglik);
/* sumloglik += loglik; */
- value2[k][i][t]=loglik;
+ value2[k][i][t]=loglik;
}
}
else
{
- LL_mod(&fam[k][i],T,vdirect[k][i],vindirect[k][(*d-m)],&theta[k][i],&nu[k][i],&loglik);
+ LL_mod(&fam[k][i],T,vdirect[k][i],vindirect[k][(*d-m)],&theta[k][i],&nu[k][i],&loglik);
/* sumloglik += loglik; */
value[k][i]=loglik;
}
- //if(cdirect[k-1][i]==1)
+ //if(cdirect[k-1][i]==1)
if(condirect[k+(*d)*i-1]==1)
{
Hfunc1(&fam[k][i],T,vdirect[k][i],vindirect[k][(*d-m)],&theta[k][i],&nu[k][i],vdirect[k-1][i]);
@@ -209,7 +209,7 @@
//if(cindirect[k-1][i]==1)
if(conindirect[k+(*d)*i-1]==1)
{
- Hfunc2(&fam[k][i],T,vindirect[k][(*d-m)],vdirect[k][i],&theta[k][i],&nu[k][i],vindirect[k-1][i]);
+ Hfunc2(&fam[k][i],T,vindirect[k][(*d-m)],vdirect[k][i],&theta[k][i],&nu[k][i],vindirect[k-1][i]);
}
}
}
@@ -218,24 +218,24 @@
sumloglik += value[k][i];
}
else
- {
- for(t=0;t<*T;t++ )
- {
+ {
+ for(t=0;t<*T;t++ )
+ {
sumsplitlog[t] += value2[k][i][t];
- }
+ }
}
- }
+ }
}
-
+
if(*seperate==0)
- {
+ {
*out = sumloglik;
}
else
{
- for(t=0;t<*T;t++ )
+ for(t=0;t<*T;t++ )
{
out[t] = sumsplitlog[t];
}
@@ -243,8 +243,8 @@
for(i=0;i<(*d);i++)
{
for(j=0;j<(*d);j++)
- {
- for(t=0;t<*T;t++ )
+ {
+ for(t=0;t<*T;t++ )
{
vv[(i+1)+(*d)*j+(*d)*(*d)*t-1]=vdirect[i][j][t];
vv2[(i+1)+(*d)*j+(*d)*(*d)*t-1]=vindirect[i][j][t];
@@ -269,18 +269,18 @@
//Free memory:
- free_matrix(x,*d);
- free_3darray(vdirect,*d,*d);
- free_matrix(theta,*d);
- free_matrix(nu,*d);
- free_intmatrix(fam,*d);
+ free_matrix(x,*d);
+ free_3darray(vdirect,*d,*d);
+ free_matrix(theta,*d);
+ free_matrix(nu,*d);
+ free_intmatrix(fam,*d);
free_matrix(value,*d);
free_3darray(value2,*d,*d);
- //free_intmatrix(mmat,*d);
- //free_intmatrix(cdirect,*d);
- //free_intmatrix(cindirect,*d);
- free_3darray(vindirect,*d,*d);
- //free_intmatrix(calc, *d);
+ //free_intmatrix(mmat,*d);
+ //free_intmatrix(cdirect,*d);
+ //free_intmatrix(cindirect,*d);
+ free_3darray(vindirect,*d,*d);
+ //free_intmatrix(calc, *d);
//free_intmatrix(mat, *d);
Free(sumsplitlog);
}
@@ -296,41 +296,41 @@
// par2 second set of parameter values (f.e. for student copulas)
// data data set for which to compute log-likelihood
// matrix an RVineMatrix in vector form
-// condirect, conindirect Matrizes which tell us where we find the right values
+// condirect, conindirect Matrizes which tell us where we find the right values
// seperate Control Parameter, do we want to seperate the likelihoods for each data point?
// calcupdate matrix which tells us for which parameters we need to redo the calculations, not newly computed values are taken from ll, vv, vv2
// Output:
// out Loglikelihood
// ll array with the contribution to LL (for each copula)
-// vv,vv2 array for the transformation operated (Hfunc)
+// vv,vv2 array for the transformation operated (Hfunc)
/////////////////////////////////////////////////////////////
/// seperate = 1 not implemented, highly experimental
-void VineLikRvine(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+void VineLikRvine(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
double* out, double* ll, double* vv, double* vv2, int* calcupdate, int* seperate)
{
int i, j, k, t, m, **fam, **cdirect, **cindirect, **mat, **mmat, **calc;
int kk=1;
double loglik=1.0, sumloglik=1.0, **x, **theta, **nu, ***vdirect, ***vindirect;
-
+
double *sumsplitlog;
sumsplitlog=(double*) Calloc(*T,double);
-
- for (t=0;t<*T;t++ )
+
+ for (t=0;t<*T;t++ )
{
sumsplitlog[t] = 0;
}
-
-
-
+
+
+
double ***value2;
value2=create_3darray(*d,*d,*T);
double **value;
value=create_matrix(*d,*d);
-
-
-
+
+
+
//Allocate memory
x = create_matrix(*d,*T);
vdirect = create_3darray(*d,*d,*T);
@@ -343,24 +343,24 @@
cindirect=create_intmatrix(*d,*d);
mat=create_intmatrix(*d,*d);
calc=create_intmatrix(*d,*d);
-
-
+
+
/* m=*d;
Rprintf("%d\n\n",m);
m=*T;
Rprintf("%d\n\n",m); */
-
+
//Initialize
k=0;
for(i=0;i<(*d);i++)
{
- for (t=0;t<*T;t++ )
+ for (t=0;t<*T;t++ )
{
x[i][t] = data[k];
k++;
}
}
-
+
k=0;
for(i=0;i<(*d);i++)
{
@@ -369,7 +369,7 @@
theta[i][j]=par[(i+1)+(*d)*j-1] ;
nu[i][j]=par2[(i+1)+(*d)*j-1] ;
mmat[i][j]=maxmat[(i+1)+(*d)*j-1] ;
- /* m=maxmat[(i+1)+(*d)*j-1];
+ /* m=maxmat[(i+1)+(*d)*j-1];
Rprintf("%d\n",m); */
mat[i][j]=matrix[(i+1)+(*d)*j-1] ;
cdirect[i][j]=condirect[(i+1)+(*d)*j-1];
@@ -378,96 +378,96 @@
if(*seperate==0){value[i][j]=ll[(i+1)+(*d)*j-1] ;}
calc[i][j]=calcupdate[(i+1)+(*d)*j-1] ;
}
- }
-
-
+ }
+
+
for(i=0;i<(*d);i++)
{
- for(j=0;j<(*d);j++){
+ for(j=0;j<(*d);j++){
for(t=0;t<*T;t++ ) {
vdirect[i][j][t]=vv[(i+1)+(*d)*j+(*d)*(*d)*t-1];
vindirect[i][j][t]=vv2[(i+1)+(*d)*j+(*d)*(*d)*t-1];
if(*seperate==1){ value2[i][j][t]=ll[(i+1)+(*d)*j+(*d)*(*d)*t-1];}
}}}
-
-
-
+
+
+
for(i=0;i<(*d);i++)
{
- for(t=0;t<*T;t++ )
+ for(t=0;t<*T;t++ )
{
vdirect[*d-1][i][t]=x[*d-1-i][t];
}
}
-
-
+
+
for(i=*d-2; i>-1; i--)
{
for(k=*d-1;k>i;k--)
- {
-
+ {
+
if(calc[k][i]==1){
m=mmat[k][i];
-
-
+
+
/* Rprintf("%d\n",m);*/
-
+
if(m == mat[k][i])
- {
+ {
if(*seperate==1){kk = 1;
for(t=0;t<*T;t++ ) {
-
- copLik_mod(&fam[k][i],&kk,&vdirect[k][i][t],&vdirect[k][(*d-m)][t],&theta[k][i],&nu[k][i],&loglik);
+
+ copLik_mod(&fam[k][i],&kk,&vdirect[k][i][t],&vdirect[k][(*d-m)][t],&theta[k][i],&nu[k][i],&loglik);
/* sumloglik += loglik; */
value2[k][i][t]=loglik;
-
+
}}else{
-
- copLik_mod(&fam[k][i],T,vdirect[k][i],vdirect[k][(*d-m)],&theta[k][i],&nu[k][i],&loglik);
+
+ copLik_mod(&fam[k][i],T,vdirect[k][i],vdirect[k][(*d-m)],&theta[k][i],&nu[k][i],&loglik);
/* sumloglik += loglik; */
value[k][i]=loglik;}
if(cdirect[k-1][i]==1) {
Hfunc1(&fam[k][i],T,vdirect[k][i],vdirect[k][*d-m],&theta[k][i],&nu[k][i],vdirect[k-1][i]);
}
if(cindirect[k-1][i]==1) {
- Hfunc2(&fam[k][i],T,vdirect[k][(*d-m)],vdirect[k][i],&theta[k][i],&nu[k][i],vindirect[k-1][i]);
+ Hfunc2(&fam[k][i],T,vdirect[k][(*d-m)],vdirect[k][i],&theta[k][i],&nu[k][i],vindirect[k-1][i]);
}
-
-
-
+
+
+
}else{
if(*seperate==1){kk = 1;
for(t=0;t<*T;t++ ) {
- copLik_mod(&fam[k][i],&kk,&vdirect[k][i][t],&vindirect[k][(*d-m)][t],&theta[k][i],&nu[k][i],&loglik);
+ copLik_mod(&fam[k][i],&kk,&vdirect[k][i][t],&vindirect[k][(*d-m)][t],&theta[k][i],&nu[k][i],&loglik);
/* sumloglik += loglik; */
value2[k][i][t]=loglik;
-
+
}}else{
-
- copLik_mod(&fam[k][i],T,vdirect[k][i],vindirect[k][(*d-m)],&theta[k][i],&nu[k][i],&loglik);
+
+ copLik_mod(&fam[k][i],T,vdirect[k][i],vindirect[k][(*d-m)],&theta[k][i],&nu[k][i],&loglik);
/* sumloglik += loglik; */
value[k][i]=loglik;}
if(cdirect[k-1][i]==1) {
Hfunc1(&fam[k][i],T,vdirect[k][i],vindirect[k][(*d-m)],&theta[k][i],&nu[k][i],vdirect[k-1][i]);
}
if(cindirect[k-1][i]==1) {
- Hfunc2(&fam[k][i],T,vindirect[k][(*d-m)],vdirect[k][i],&theta[k][i],&nu[k][i],vindirect[k-1][i]);
+ Hfunc2(&fam[k][i],T,vindirect[k][(*d-m)],vdirect[k][i],&theta[k][i],&nu[k][i],vindirect[k-1][i]);
}
-
+
}
}
if(*seperate==0){sumloglik = sumloglik*value[k][i];}else{ for(t=0;t<*T;t++ ) { sumsplitlog[t] += value2[k][i][t];} };
- }
+ }
}
-
-
-
-
-
+
+
+
+
+
if(*seperate==0){ *out = sumloglik;}else{for(t=0;t<*T;t++ ) {out[t] = sumsplitlog[t];}};
for(i=0;i<(*d);i++)
{
- for(j=0;j<(*d);j++){
+ for(j=0;j<(*d);j++){
for(t=0;t<*T;t++ ) {
vv[(i+1)+(*d)*j+(*d)*(*d)*t-1]=vdirect[i][j][t];
vv2[(i+1)+(*d)*j+(*d)*(*d)*t-1]=vindirect[i][j][t];
@@ -478,8 +478,8 @@
for(j=0;j<(*d);j++)
{
if(*seperate==0){ll[(i+1)+(*d)*j-1]=value[i][j];}}}
-
-
+
+
//Free memory:
free_matrix(x,*d); free_3darray(vdirect,*d,*d); free_matrix(theta,*d); free_matrix(nu,*d); free_intmatrix(fam,*d); free_matrix(value,*d);free_3darray(value2,*d,*d);
free_intmatrix(mmat,*d); free_intmatrix(cdirect,*d); free_intmatrix(cindirect,*d); free_3darray(vindirect,*d,*d); free_intmatrix(calc, *d); free_intmatrix(mat, *d);
@@ -494,11 +494,11 @@
// PCC for R-vine
////////////////////////////////////////////
-void SimulateRVine(int* T, int* d, int* family, int* maxmat, int* matrix, int* conindirect, double* par, double* par2, double* out)
+void SimulateRVine(int* T, int* d, int* family, int* maxmat, int* matrix, int* conindirect, double* par, double* par2, double* out, double* U, int* takeU)
{
int i, j, k, m, **fam, **cindirect, **mat, **mmat, **fam2, **cindirect2, **mat2, **mmat2;
- double **theta, **nu, **theta2, **nu2, ***vdirect, ***vindirect;
-
+ double **theta, **nu, **theta2, **nu2, ***vdirect, ***vindirect, **U2;
+
//Allocate memory
theta=create_matrix(*d,*d);
nu=create_matrix(*d,*d);
@@ -514,6 +514,7 @@
mat2=create_intmatrix(*d,*d);
vdirect = create_3darray(*d,*d,*T);
vindirect = create_3darray(*d,*d,*T);
+ U2 = create_matrix(*T, *d);
//Initialize random number generator:
GetRNGstate();
@@ -521,18 +522,19 @@
//Initialize
k=0;
for(i=0;i<(*d);i++)
- {
- for(j=0;j<(*d);j++)
- {
- theta2[i][j]=par[(i+1)+(*d)*j-1] ;
- nu2[i][j]=par2[(i+1)+(*d)*j-1] ;
- mmat2[i][j]=maxmat[(i+1)+(*d)*j-1] ;
- mat2[i][j]=matrix[(i+1)+(*d)*j-1] ;
- cindirect2[i][j]=conindirect[(i+1)+(*d)*j-1] ;
- fam2[i][j]=family[(i+1)+(*d)*j-1] ;
+ {
+ for(j=0;j<(*d);j++)
+ {
+ theta2[i][j]=par[(i+1)+(*d)*j-1] ;
+ nu2[i][j]=par2[(i+1)+(*d)*j-1] ;
+ mmat2[i][j]=maxmat[(i+1)+(*d)*j-1] ;
+ mat2[i][j]=matrix[(i+1)+(*d)*j-1] ;
+ cindirect2[i][j]=conindirect[(i+1)+(*d)*j-1] ;
+ fam2[i][j]=family[(i+1)+(*d)*j-1] ;
}
- }
-
+ }
+ for(j=0;j<(*d);j++) for(i=0;i<(*T);i++) U2[i][j]=U[(*T)*j+i]; // (T [=N], d)-matrix
+
// Matrizen rotieren für den Algo
for(i=0;i<(*d);i++)
{
@@ -548,12 +550,12 @@
}
free_matrix(theta2,*d);
- free_matrix(nu2,*d);
+ free_matrix(nu2,*d);
free_intmatrix(fam2,*d);
- free_intmatrix(mmat2,*d);
- free_intmatrix(cindirect2,*d);
+ free_intmatrix(mmat2,*d);
+ free_intmatrix(cindirect2,*d);
free_intmatrix(mat2, *d);
-
+
/*
Declare variable to hold seconds on clock.
*/
@@ -571,9 +573,10 @@
// Der eigentliche Algo
- for(j=0;j<*T;j++)
- {
- for(i=0;i<*d;i++) vdirect[i][i][j] = runif(0,1);
+ for(j=0;j<*T;j++) // sample size
+ {
+ if(*takeU == 1) for(i=0;i<*d;i++) vdirect[i][i][j] = U2[j][i]; // j = 'sample size'; i = 'copula dimension'
+ else for(i=0;i<*d;i++) vdirect[i][i][j] = runif(0,1);
vindirect[0][0][j] = vdirect[0][0][j];
}
@@ -607,8 +610,8 @@
}
}
}
-
+
k=0;
for(i=0;i<(*d);i++)
{
@@ -621,12 +624,13 @@
//Free memory:
free_matrix(theta,*d);
- free_matrix(nu,*d);
+ free_matrix(nu,*d);
free_intmatrix(fam,*d);
- free_intmatrix(mmat,*d);
- free_intmatrix(cindirect,*d);
+ free_intmatrix(mmat,*d);
+ free_intmatrix(cindirect,*d);
free_intmatrix(mat, *d);
free_3darray(vdirect,*d,*d);
free_3darray(vindirect,*d,*d);
+ free_matrix(U2,*T);
PutRNGstate();
}
Mehr Informationen über die Mailingliste Vinecopula-commits