[Vinecopula-commits] r89 - pkg/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
So Mär 29 18:07:38 CEST 2015
Author: ulf
Date: 2015-03-29 18:07:38 +0200 (Sun, 29 Mar 2015)
New Revision: 89
Modified:
pkg/src/cdvine.c
pkg/src/deriv.c
Log:
I started to comment my C-functions
Modified: pkg/src/cdvine.c
===================================================================
--- pkg/src/cdvine.c 2015-03-28 13:29:18 UTC (rev 88)
+++ pkg/src/cdvine.c 2015-03-29 16:07:38 UTC (rev 89)
@@ -10,12 +10,14 @@
**
*/
-#include "include/vine.h"
-#include "include/memoryhandling.h"
-#include "include/likelihood.h"
-#include "include/cdvine.h"
-#include "include/hfunc.h"
+// 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
@@ -24,13 +26,17 @@
//////////////////////////////////////////////////////////////
-// Function to simulate from a pair-copula construction (vine)
+// 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 (1=gaussian, 2=student, 3=clayton, 4=gumbel, 5=frank, 6=joe, 7=BB1)
+// 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)
@@ -38,7 +44,7 @@
int i, j, in=1, k, **fam;
double *w, **v, t, **theta, **x, **ny;
- GetRNGstate();
+ GetRNGstate(); //Init random number generator
//Allocate memory:
w = Calloc((*d+1),double);
@@ -48,6 +54,9 @@
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++)
{
@@ -59,14 +68,14 @@
k ++;
}
}
- //Simulate:
+ //Simulate: (it follows the theoretical algorithm)
if(*type==1) //Canonical vine
{
- for(j=1;j<=*n;j++)
+ 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++)
+ for(i=2;i<=*d;i++) // run over all dimensions (cols)
{
t = w[i];
for(k=i-1;k>=1;k--)
@@ -134,7 +143,7 @@
k ++;
}
}
- PutRNGstate();
+ 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);
}
@@ -143,16 +152,16 @@
//////////////////////////////////////////////////////////////
-// Function to compute -log-likelihood for the pair-copula construction (vine)
+// 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: only student // (1=gaussian, 2=student, 3=clayton, 4=gumbel, 5=frank, 7=BB1, 8=BB7)
-// par parameter values (at least d*(d-1)/2 parameters
+// 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 Loglikelihood
+// out Log-likelihood
// ll array with the contribution to LL (for each copula)
// vv array for the transformation operated (Hfunc)
/////////////////////////////////////////////////////////////
@@ -184,7 +193,7 @@
{
for (t=0;t<=*T-1;t++ )
{
- x[i][t] = data[k];
+ x[i][t] = data[k]; //transform the data back into a 2-dim array
k++;
}
}
@@ -195,7 +204,7 @@
{
theta[i][j] = par[k];
fam[i][j] = family[k];
- nu[i][j] = par[*d*(*d-1)/2+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++;
}
}
@@ -207,9 +216,10 @@
//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);
- sumloglik += loglik;
- ll[kk] = loglik;
+ 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)
{
@@ -297,7 +307,7 @@
for(i=1;i<=(*d-k);i++)
for(t=0;t<*T;t++)
{
- vv[kk] = v[k][i][t];
+ vv[kk] = v[k][i][t]; // transformation from a 3-dim array to a vector
++kk;
}
}
@@ -324,18 +334,19 @@
}
//////////////////////////////////////////////////////////////
-// Function to compute -log-likelihood for the pair-copula construction (vine)
+// 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: only student // (1=gaussian, 2=student, 3=clayton, 4=gumbel, 5=frank)
+// 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)
/////////////////////////////////////////////////////////////
Modified: pkg/src/deriv.c
===================================================================
--- pkg/src/deriv.c 2015-03-28 13:29:18 UTC (rev 88)
+++ pkg/src/deriv.c 2015-03-29 16:07:38 UTC (rev 89)
@@ -23,11 +23,21 @@
/////////////////////////////////////////////////////////////
//
// Ableitung der Copula nach dem Parameter
+// Derivative of bivariate copulas with respect to the (first) parameter
//
+// Input:
+// u,v copula arguments (data vectors)
+// n length of u,v
+// param parameter vector (par,par2)
+// copula copula family
+//
+// Output:
+// out derivative
/////////////////////////////////////////////////////////////
void diffPDF_mod(double* u, double* v, int* n, double* param, int* copula, double* out)
{
+ // for the rotated copulas we need some help variables
double* negv;
double* negu;
double* nparam;
@@ -39,12 +49,12 @@
nparam[1]=-param[1];
int i;
-if((*copula)==43)
+if((*copula)==43) // special copula; all rotations of Clayton are combined in one copula; the correct rotation is derived automatically
{
ncopula=3;
if(param[0] > 0){
nparam[0]=2*(param[0])/(1-param[0]);
- diffPDF(u, v, n, nparam, &ncopula, out);
+ diffPDF(u, v, n, nparam, &ncopula, out); // derivative function
for(i=0;i<*n;i++){out[i]=out[i]*2/pow(1-param[0],2);}
}else{
nparam[0]=-2*(param[0])/(1+param[0]);
@@ -66,22 +76,22 @@
for(i=0;i<*n;i++){out[i]=-out[i]/pow(1+param[0],2); }
}
}else{
-
- if(((*copula==23) | (*copula==24) | (*copula==26) | (*copula==27) | (*copula==28) | (*copula==29) | (*copula==30))) // 90? rotated copulas
+// for the rotation see the master thesis of Jakob Stoeber
+ if(((*copula==23) | (*copula==24) | (*copula==26) | (*copula==27) | (*copula==28) | (*copula==29) | (*copula==30))) // 90 rotated copulas
{
ncopula = (*copula)-20;
for (i = 0; i < *n; ++i) {negv[i] = 1 - v[i];}
diffPDF(u, negv, n, nparam, &ncopula, out);
for(i=0;i<*n;i++){out[i]=-out[i];}
}
- else if(((*copula==33) | (*copula==34) | (*copula==36) | (*copula==37) | (*copula==38) | (*copula==39) | (*copula==40))) // 270? rotated copulas
+ else if(((*copula==33) | (*copula==34) | (*copula==36) | (*copula==37) | (*copula==38) | (*copula==39) | (*copula==40))) // 270 rotated copulas
{
ncopula = (*copula)-30;
for (i = 0; i < *n; ++i) {negu[i] = 1 - u[i];}
diffPDF(negu, v, n, nparam, &ncopula, out);
for(i=0;i<*n;i++){out[i]=-out[i];}
}
- else if(((*copula==13) | (*copula==14) | (*copula==16) | (*copula==17) | (*copula==18) | (*copula==19) | (*copula==20))) // 180? rotated copulas
+ else if(((*copula==13) | (*copula==14) | (*copula==16) | (*copula==17) | (*copula==18) | (*copula==19) | (*copula==20))) // 180 rotated copulas
{
ncopula = (*copula)-10;
for (i = 0; i < *n; ++i)
@@ -93,7 +103,7 @@
}
else
{
- diffPDF(u, v, n, param, copula, out);
+ diffPDF(u, v, n, param, copula, out); // eigentliche Ableitungsfunktion
}
}
free(negv);
@@ -104,6 +114,21 @@
+//////////////////////////////////////////////////
+// Derivative of bivariate copulas with respect to the parameter (standard copula form without rotations, see above)
+//
+// Input:
+// u,v copula arguments (data vectors)
+// n length of u,v
+// param parameter vector (par,par2)
+// copula copula family (1,3,4,5,6)
+//
+// Output:
+// out derivative
+//
+// Reference: Schepsmeier and Stoeber (2012, 2013)
+/////////////////////////////////////////////////////////////
+
void diffPDF(double* u, double* v, int* n, double* param, int* copula, double* out)
{
int j;
@@ -116,11 +141,11 @@
for(j=0;j<*n;j++)
{
- if(*copula==0)
+ if(*copula==0) // independence copulas
{
out[j]=0;
}
- else if(*copula==1)
+ else if(*copula==1) // gauss
{
t1 = qnorm(u[j],0.0,1.0,1,0);
t2 = qnorm(v[j],0.0,1.0,1,0);
@@ -135,6 +160,7 @@
t24 = sqrt(t8);
out[j] = (-2.0*(theta*t3-t1*t2)*t9-t15/(t8*t8)*theta)*t22/t24+t22/t24/t8*theta;
}
+ // t-copula is separate; very complicated
else if(*copula==3)
{
t1 = u[j]*v[j];
@@ -240,7 +266,16 @@
////////////////////////////////////////////////////////////////////
//
// 1. Ableitung von c nach u
+// First derivative of the bivariate copula density with respect to u (first argument)
+// Input:
+// u,v copula arguments (data vectors)
+// n length of u,v
+// param parameter vector (par,par2)
+// copula copula family
//
+// Output:
+// out derivative
+//
////////////////////////////////////////////////////////////////////
void diffPDF_u_mod(double* u, double* v, int* n, double* param, int* copula, double* out)
@@ -316,6 +351,22 @@
free(nparam);
}
+
+////////////////////////////////////////////////////////////////////
+//
+// 1. Ableitung von c nach u (eigentliche Funktion ohne die Rotationen)
+// First derivative of the bivariate copula density with respect to u (first argument)
+// Input:
+// u,v copula arguments (data vectors)
+// n length of u,v
+// param parameter vector (par,par2)
+// copula copula family (1,2,3,4,5,6)
+//
+// Output:
+// out derivative
+//
+////////////////////////////////////////////////////////////////////
+
void diffPDF_u(double* u, double* v, int* n, double* param, int* copula, double* out)
{
int j, k=1;
@@ -350,7 +401,7 @@
}
else if(*copula==2)
{
- diffPDF_u_tCopula_new(&u[j], &v[j], &k, param, copula, &out[j]);
+ diffPDF_u_tCopula_new(&u[j], &v[j], &k, param, copula, &out[j]); // special function for t-copula
}
else if(*copula==3)
{
@@ -436,7 +487,16 @@
////////////////////////////////////////////////////////////////////
//
// 1. Ableitung von c nach v
+// First derivative of the bivariate copula density with respect to v (second argument)
+// Input:
+// u,v copula arguments (data vectors)
+// n length of u,v
+// param parameter vector (par,par2)
+// copula copula family (1,2,3,4,5,6)
//
+// Output:
+// out derivative
+//
////////////////////////////////////////////////////////////////////
void diffPDF_v_mod(double* u, double* v, int* n, double* param, int* copula, double* out)
@@ -480,7 +540,7 @@
{
ncopula = (*copula)-20;
for (i = 0; i < *n; ++i) {negv[i] = 1 - v[i];}
- diffPDF_u(negv, u, n, nparam, &ncopula, out);
+ diffPDF_u(negv, u, n, nparam, &ncopula, out); // we can use again the function for the derivative of c wrt u but change the arguments
for(i=0;i<*n;i++){out[i]=-out[i];}
}
else if(((*copula==33) | (*copula==34) | (*copula==36) | (*copula==37) | (*copula==38) | (*copula==39) | (*copula==40))) // 270? rotated copulas
Mehr Informationen über die Mailingliste Vinecopula-commits