[Vinecopula-commits] r136 - in pkg: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Do Sep 17 18:50:01 CEST 2015
Author: tnagler
Date: 2015-09-17 18:50:00 +0200 (Thu, 17 Sep 2015)
New Revision: 136
Modified:
pkg/R/BiCopTau2Par.r
pkg/src/likelihood.c
Log:
* allow for zero parameter for Frank and Clayton BiCopPDF (so far only when check.pars = FALSE) and BiCopTau2Par
Modified: pkg/R/BiCopTau2Par.r
===================================================================
--- pkg/R/BiCopTau2Par.r 2015-09-17 12:24:05 UTC (rev 135)
+++ pkg/R/BiCopTau2Par.r 2015-09-17 16:50:00 UTC (rev 136)
@@ -34,32 +34,30 @@
} else if (family %in% 1:2) {
par <- sin(pi * tau/2)
} else if (family %in% c(3, 13)) {
- if (tau <= 0)
- stop("Clayton copula cannot be used for tau<=0.")
+ if (tau < 0)
+ stop("Clayton copula cannot be used for tau<0.")
par <- 2 * tau/(1 - tau)
} else if (family %in% c(4, 14)) {
if (tau < 0)
stop("Gumbel copula cannot be used for tau<0.")
par <- 1/(1 - tau)
} else if (family == 5) {
- if (tau == 0)
- stop("Frank copula cannot be used for tau=0.")
- par <- Frank.itau.JJ(tau)
+ par <- if (tau == 0) 0 else Frank.itau.JJ(tau)
} else if (family %in% c(6, 16)) {
- if (tau <= 0)
- stop("Joe copula cannot be used for tau<=0.")
+ if (tau < 0)
+ stop("Joe copula cannot be used for tau<0.")
par <- Joe.itau.JJ(tau)
} else if (family %in% c(23, 33)) {
- if (tau >= 0)
- stop("Rotated Clayton copula cannot be used for tau>=0.")
+ if (tau > 0)
+ stop("Rotated Clayton copula cannot be used for tau>0.")
par <- 2 * tau/(1 + tau)
} else if (family %in% c(24, 34)) {
if (tau > 0)
stop("Rotated Gumbel copula cannot be used for tau>0.")
par <- -(1/(1 + tau))
} else if (family %in% c(26, 36)) {
- if (tau >= 0)
- stop("Rotated Joe copula cannot be used for tau>=0.")
+ if (tau > 0)
+ stop("Rotated Joe copula cannot be used for tau>0.")
par <- -Joe.itau.JJ(-tau)
} else if (family %in% c(41, 51)) {
par <- ipsA.tau2cpar(tau)
Modified: pkg/src/likelihood.c
===================================================================
--- pkg/src/likelihood.c 2015-09-17 12:24:05 UTC (rev 135)
+++ pkg/src/likelihood.c 2015-09-17 16:50:00 UTC (rev 136)
@@ -1,14 +1,14 @@
/*
-** likelihood.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.
-**
-*/
+ ** likelihood.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 "include/vine.h"
#include "include/memoryhandling.h"
@@ -35,57 +35,57 @@
// out outout
//////////////////////////////////////////////////////////////
/*
-void gen(double* u, int* n, double* param, int* copula, double* out)
-{
- int j;
- double *h;
- h = Calloc(*n,double);
+ void gen(double* u, int* n, double* param, int* copula, double* out)
+ {
+ int j;
+ double *h;
+ h = Calloc(*n,double);
+
+ for(j=0;j<*n;j++)
+ {
+ if(u[j]==0) h[j] = 0;
+ else if (u[j]==1) h[j] = u[j];
+ else
+ {
+ if(*copula==3) //Clayton
+ {
+ h[j] = 1/param[0]*(pow(u[j],(-param[0]))-1);
+ }
+ if(*copula==4) //Gumbel
+ {
+ h[j] = pow((-log(u[j])),param[0]);
+ }
+ if(*copula==5) //Frank
+ {
+ h[j] = -log((exp(-param[0]*u[j])-1)/(exp(-param[0])-1));
+ }
+ if(*copula==6) //Joe
+ {
+ h[j] = -log(1-pow((1-u[j]),param[0]));
+ }
+ if(*copula==7) //BB1
+ {
+ h[j] = pow((pow(u[j],(-param[0]))-1),param[1]);
+ }
+ else if(*copula==8) //BB6
+ {
+ h[j] = pow((-log(-pow(1-u[j],param[0])+1)),param[1]);
+ }
+ else if(*copula==9) //BB7
+ {
+ h[j] = pow(1-pow(1-u[j],param[0]),-param[1])-1;
+ }
+ else if(*copula==10) //BB8
+ {
+ h[j] = -log( (1-pow(1-param[1]*u[j],param[0])) / (1-pow(1-param[1],param[0])) );
+ }
+ }
+ out[j]=h[j];
+ }
+ Free(h);
+ }
+ */
- for(j=0;j<*n;j++)
- {
- if(u[j]==0) h[j] = 0;
- else if (u[j]==1) h[j] = u[j];
- else
- {
- if(*copula==3) //Clayton
- {
- h[j] = 1/param[0]*(pow(u[j],(-param[0]))-1);
- }
- if(*copula==4) //Gumbel
- {
- h[j] = pow((-log(u[j])),param[0]);
- }
- if(*copula==5) //Frank
- {
- h[j] = -log((exp(-param[0]*u[j])-1)/(exp(-param[0])-1));
- }
- if(*copula==6) //Joe
- {
- h[j] = -log(1-pow((1-u[j]),param[0]));
- }
- if(*copula==7) //BB1
- {
- h[j] = pow((pow(u[j],(-param[0]))-1),param[1]);
- }
- else if(*copula==8) //BB6
- {
- h[j] = pow((-log(-pow(1-u[j],param[0])+1)),param[1]);
- }
- else if(*copula==9) //BB7
- {
- h[j] = pow(1-pow(1-u[j],param[0]),-param[1])-1;
- }
- else if(*copula==10) //BB8
- {
- h[j] = -log( (1-pow(1-param[1]*u[j],param[0])) / (1-pow(1-param[1],param[0])) );
- }
- }
- out[j]=h[j];
- }
- Free(h);
-}
-*/
-
//////////////////////////////////////////////////////////////
// Inverse generator of BB1, BB6, BB7 and BB8
// Input:
@@ -96,57 +96,57 @@
// out outout
//////////////////////////////////////////////////////////////
/*
-void genInv(double* u, int* n, double* param, int* copula, double* out)
-{
- int j;
- double *h;
- h = Calloc(*n,double);
+ void genInv(double* u, int* n, double* param, int* copula, double* out)
+ {
+ int j;
+ double *h;
+ h = Calloc(*n,double);
+
+ for(j=0;j<*n;j++)
+ {
+ if(u[j]==0) h[j] = 0;
+ else if (u[j]==1) h[j] = u[j];
+ else
+ {
+ if(*copula==3) //Clayton
+ {
+ h[j] = pow((1+param[0]*u[j]),(-1/param[0]));
+ }
+ if(*copula==4) //Gumbel
+ {
+ h[j] = exp(-pow(u[j],1/param[0]));
+ }
+ if(*copula==5) //Frank
+ {
+ h[j] = -1/param[0]*log(1-exp(-u[j])*(1-exp(-param[0])));
+ }
+ if(*copula==6) //Joe
+ {
+ h[j] = 1-pow((1-exp(-u[j])),1/param[0]);
+ }
+ if(*copula==7) //BB1
+ {
+ h[j] = pow(1+pow(u[j],1/param[1]),(-1/param[0]));
+ }
+ else if(*copula==8) //BB6
+ {
+ h[j] = 1-pow(1-exp(-pow(u[j],1/param[1])),1/param[0]);
+ }
+ else if(*copula==9) //BB7
+ {
+ h[j] = 1-pow(1-pow(1+u[j],-1/param[1]),(1/param[0]));
+ }
+ else if(*copula==10) //BB8
+ {
+ h[j] = 1/param[1] * ( 1-pow(1-(1-pow(1-param[1],param[0]))*exp(-u[j]),1/param[0]) );
+ }
+ }
+ out[j]=h[j];
+ }
+ Free(h);
+ }
+ */
- for(j=0;j<*n;j++)
- {
- if(u[j]==0) h[j] = 0;
- else if (u[j]==1) h[j] = u[j];
- else
- {
- if(*copula==3) //Clayton
- {
- h[j] = pow((1+param[0]*u[j]),(-1/param[0]));
- }
- if(*copula==4) //Gumbel
- {
- h[j] = exp(-pow(u[j],1/param[0]));
- }
- if(*copula==5) //Frank
- {
- h[j] = -1/param[0]*log(1-exp(-u[j])*(1-exp(-param[0])));
- }
- if(*copula==6) //Joe
- {
- h[j] = 1-pow((1-exp(-u[j])),1/param[0]);
- }
- if(*copula==7) //BB1
- {
- h[j] = pow(1+pow(u[j],1/param[1]),(-1/param[0]));
- }
- else if(*copula==8) //BB6
- {
- h[j] = 1-pow(1-exp(-pow(u[j],1/param[1])),1/param[0]);
- }
- else if(*copula==9) //BB7
- {
- h[j] = 1-pow(1-pow(1+u[j],-1/param[1]),(1/param[0]));
- }
- else if(*copula==10) //BB8
- {
- h[j] = 1/param[1] * ( 1-pow(1-(1-pow(1-param[1],param[0]))*exp(-u[j]),1/param[0]) );
- }
- }
- out[j]=h[j];
- }
- Free(h);
-}
-*/
-
//////////////////////////////////////////////////////////////
// First derivative of the generator of BB1, BB2 and BB7
// Input:
@@ -157,40 +157,40 @@
// out outout
//////////////////////////////////////////////////////////////
/*
-void genDrv(double* u, int* n, double* param, int* copula, double* out)
-{
- int j;
- double *h;
- h = Calloc(*n,double);
-
- for(j=0;j<*n;j++)
- {
- if(u[j]==0) h[j] = 0;
- else if (u[j]==1) h[j] = u[j];
- else
- {
- if(*copula==7) //BB1
- {
- h[j] = -(param[0]*param[1])*pow(pow(u[j],-param[0])-1,param[1]-1)*pow(u[j],-1-param[0]);
- }
- else if(*copula==8) //BB6
- {
- h[j] = pow(-log(-pow(1-u[j],param[0])+1),param[1]-1) * param[1] * pow(1-u[j],param[0]-1) * param[0] / (pow(1-u[j],param[0])-1);
- }
- else if(*copula==9) //BB7
- {
- h[j] = -(param[0]*param[1])*pow(1-u[j],param[0]-1)*pow(1-pow(1-u[j],param[0]),-1-param[1]);
- }
- else if(*copula==10) //BB8
- {
- h[j] = -param[0] * param[1] * ( pow(1-param[1]*u[j],param[0]-1) ) / ( 1-pow(1-param[1]*u[j],param[0]) );
- }
- }
- out[j]=h[j];
- }
- Free(h);
-}
-*/
+ void genDrv(double* u, int* n, double* param, int* copula, double* out)
+ {
+ int j;
+ double *h;
+ h = Calloc(*n,double);
+
+ for(j=0;j<*n;j++)
+ {
+ if(u[j]==0) h[j] = 0;
+ else if (u[j]==1) h[j] = u[j];
+ else
+ {
+ if(*copula==7) //BB1
+ {
+ h[j] = -(param[0]*param[1])*pow(pow(u[j],-param[0])-1,param[1]-1)*pow(u[j],-1-param[0]);
+ }
+ else if(*copula==8) //BB6
+ {
+ h[j] = pow(-log(-pow(1-u[j],param[0])+1),param[1]-1) * param[1] * pow(1-u[j],param[0]-1) * param[0] / (pow(1-u[j],param[0])-1);
+ }
+ else if(*copula==9) //BB7
+ {
+ h[j] = -(param[0]*param[1])*pow(1-u[j],param[0]-1)*pow(1-pow(1-u[j],param[0]),-1-param[1]);
+ }
+ else if(*copula==10) //BB8
+ {
+ h[j] = -param[0] * param[1] * ( pow(1-param[1]*u[j],param[0]-1) ) / ( 1-pow(1-param[1]*u[j],param[0]) );
+ }
+ }
+ out[j]=h[j];
+ }
+ Free(h);
+ }
+ */
//////////////////////////////////////////////////////////////
// Second derivative of the generator of BB1, BB6, BB7 and BB8
@@ -202,47 +202,47 @@
// out outout
//////////////////////////////////////////////////////////////
/*
-void genDrv2(double* u, int* n, double* param, int* copula, double* out)
-{
- int j;
- double *h;
- h = Calloc(*n,double);
-
- for(j=0;j<*n;j++)
- {
- if(u[j]==0) h[j] = 0;
- else if (u[j]==1) h[j] = u[j];
- else
- {
- if(*copula==7) //BB1
- {
- h[j] = param[0]*param[1]*pow(u[j],-2-param[0])*pow(pow(u[j],-param[0])-1,param[1]-2)*((1+param[0]*param[1])*pow(u[j],-param[0])-param[0]-1);
- }
- else if(*copula==8) //BB6
- {
- h[j] = ( param[0]*param[1] * (
- pow(-log(-pow(1-u[j],param[0])+1),param[1]-2) * pow(1-u[j],2*param[0]-2) * param[0]*param[1] -
- pow(-log(-pow(1-u[j],param[0])+1),param[1]-2) * pow(1-u[j],2*param[0]-2) * param[0] -
- pow(-log(-pow(1-u[j],param[0])+1),param[1]-1) * pow(1-u[j],param[0]-2) * param[0] -
- pow(-log(-pow(1-u[j],param[0])+1),param[1]-1) * pow(1-u[j],2*param[0]-2) +
- pow(-log(-pow(1-u[j],param[0])+1),param[1]-1) * pow(1-u[j],param[0]-2) ) ) / pow(pow(1-u[j],param[0])-1,2);
- }
- else if(*copula==9) //BB7
- {
- h[j] = param[0]*param[1]*pow(1-u[j],param[0]-2)*pow(1-pow(1-u[j],param[0]),-2-param[1])*((1+param[0]*param[1])*pow(1-u[j],param[0])+param[0]-1);
- }
- else if(*copula==10) //BB8
- {
- h[j] = ( pow(param[1],2) * param[0] *
- (pow(1-u[j]*param[1],param[0]-2) * param[0] + pow(1-u[j]*param[1],2*param[0]-2) - pow(1-u[j]*param[1],param[0]-2)) ) /
- ( pow(pow(1-u[j]*param[1],param[0])-1,2) );
- }
- }
- out[j]=h[j];
- }
- Free(h);
-}
-*/
+ void genDrv2(double* u, int* n, double* param, int* copula, double* out)
+ {
+ int j;
+ double *h;
+ h = Calloc(*n,double);
+
+ for(j=0;j<*n;j++)
+ {
+ if(u[j]==0) h[j] = 0;
+ else if (u[j]==1) h[j] = u[j];
+ else
+ {
+ if(*copula==7) //BB1
+ {
+ h[j] = param[0]*param[1]*pow(u[j],-2-param[0])*pow(pow(u[j],-param[0])-1,param[1]-2)*((1+param[0]*param[1])*pow(u[j],-param[0])-param[0]-1);
+ }
+ else if(*copula==8) //BB6
+ {
+ h[j] = ( param[0]*param[1] * (
+ pow(-log(-pow(1-u[j],param[0])+1),param[1]-2) * pow(1-u[j],2*param[0]-2) * param[0]*param[1] -
+ pow(-log(-pow(1-u[j],param[0])+1),param[1]-2) * pow(1-u[j],2*param[0]-2) * param[0] -
+ pow(-log(-pow(1-u[j],param[0])+1),param[1]-1) * pow(1-u[j],param[0]-2) * param[0] -
+ pow(-log(-pow(1-u[j],param[0])+1),param[1]-1) * pow(1-u[j],2*param[0]-2) +
+ pow(-log(-pow(1-u[j],param[0])+1),param[1]-1) * pow(1-u[j],param[0]-2) ) ) / pow(pow(1-u[j],param[0])-1,2);
+ }
+ else if(*copula==9) //BB7
+ {
+ h[j] = param[0]*param[1]*pow(1-u[j],param[0]-2)*pow(1-pow(1-u[j],param[0]),-2-param[1])*((1+param[0]*param[1])*pow(1-u[j],param[0])+param[0]-1);
+ }
+ else if(*copula==10) //BB8
+ {
+ h[j] = ( pow(param[1],2) * param[0] *
+ (pow(1-u[j]*param[1],param[0]-2) * param[0] + pow(1-u[j]*param[1],2*param[0]-2) - pow(1-u[j]*param[1],param[0]-2)) ) /
+ ( pow(pow(1-u[j]*param[1],param[0])-1,2) );
+ }
+ }
+ out[j]=h[j];
+ }
+ Free(h);
+ }
+ */
//////////////////////////////////////////////////////////////
// Copula of BB1, BB6, BB7 and BB8
@@ -255,27 +255,27 @@
// out outout
//////////////////////////////////////////////////////////////
/*
-void copCdf(double* u, double* v, int* n, double* param, int* copula, double* out)
-{
- int j;
- double *out1;
- double *out2;
- double *out3;
- out1 = Calloc(*n,double);
- out2 = Calloc(*n,double);
- out3 = Calloc(*n,double);
- gen(u, n, param, copula, out1);
- gen(v, n, param, copula, out2);
- for(j=0;j<*n;j++)
- {
- out3[j]=out1[j]+out2[j];
- }
- genInv(out3 , n, param, copula, out);
- Free(out1);
- Free(out2);
- Free(out3);
-}
-*/
+ void copCdf(double* u, double* v, int* n, double* param, int* copula, double* out)
+ {
+ int j;
+ double *out1;
+ double *out2;
+ double *out3;
+ out1 = Calloc(*n,double);
+ out2 = Calloc(*n,double);
+ out3 = Calloc(*n,double);
+ gen(u, n, param, copula, out1);
+ gen(v, n, param, copula, out2);
+ for(j=0;j<*n;j++)
+ {
+ out3[j]=out1[j]+out2[j];
+ }
+ genInv(out3 , n, param, copula, out);
+ Free(out1);
+ Free(out2);
+ Free(out3);
+ }
+ */
//////////////////////////////////////////////////////////////
// Copula density of BB1, BB6, BB7 and BB8
@@ -288,183 +288,183 @@
// out outout
//////////////////////////////////////////////////////////////
/*
-void copPdf(double* u, double* v, int* n, double* param, int* copula, double* out)
-{
- int j;
- double *out1, *out2, *out3, *out4, *out5;
- out1 = Calloc(*n,double);
- out2 = Calloc(*n,double);
- out3 = Calloc(*n,double);
- out4 = Calloc(*n,double);
- out5 = Calloc(*n,double);
- copCdf(u,v,n,param,copula,out1);
- genDrv2(out1,n,param,copula,out2);
- genDrv(u,n,param,copula,out3);
- genDrv(v,n,param,copula,out4);
- genDrv(out1,n,param,copula,out5);
- for(j=0;j<*n;j++)
- {
- out[j]=-(out2[j]*out3[j]*out4[j])/pow(out5[j],3);
- }
- Free(out1);
- Free(out2);
- Free(out3);
- Free(out4);
- Free(out5);
-}
-*/
+ void copPdf(double* u, double* v, int* n, double* param, int* copula, double* out)
+ {
+ int j;
+ double *out1, *out2, *out3, *out4, *out5;
+ out1 = Calloc(*n,double);
+ out2 = Calloc(*n,double);
+ out3 = Calloc(*n,double);
+ out4 = Calloc(*n,double);
+ out5 = Calloc(*n,double);
+ copCdf(u,v,n,param,copula,out1);
+ genDrv2(out1,n,param,copula,out2);
+ genDrv(u,n,param,copula,out3);
+ genDrv(v,n,param,copula,out4);
+ genDrv(out1,n,param,copula,out5);
+ for(j=0;j<*n;j++)
+ {
+ out[j]=-(out2[j]*out3[j]*out4[j])/pow(out5[j],3);
+ }
+ Free(out1);
+ Free(out2);
+ Free(out3);
+ Free(out4);
+ Free(out5);
+ }
+ */
///////////////////////////////////////////////////////
// New
double log1mexp(double a)
{
- double result;
-if (a<log(2)) {
- result=log(-expm1(-a));
-}else{
- result=log1p(-exp(-a));
-}
- return result;
+ double result;
+ if (a<log(2)) {
+ result=log(-expm1(-a));
+ }else{
+ result=log1p(-exp(-a));
+ }
+ return result;
}
void archCDF(double* u, double* v, int* n, double* param, int* copula, double* out)
{
- int j;
- double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14;
-
- for(j=0;j<*n;j++)
- {
- if(u[j]>UMAX && v[j]>UMAX){ out[j]=1;}
- else if(u[j]>UMAX){ out[j]=v[j];}
- else if(v[j]>UMAX){ out[j]=u[j];}
- else if(u[j]<UMIN || v[j]<UMIN){ out[j]=0;}
- else
- {
- if(*copula==3) //Clayton
- {
- t1 = pow(u[j],-param[0]);
- t2 = pow(v[j],-param[0]);
- t3 = t1+t2-1;
- out[j] = pow(t3,-1/param[0]);
- }
- else if(*copula==4) //Gumbel
- {
- t1 = -log(u[j]);
- t2 = -log(v[j]);
- t3 = pow(t1,param[0]);
- t4 = pow(t2,param[0]);
- t5 = t3+t4;
- t6 = -pow(t5,1/param[0]);
- out[j] = exp(t6);
- }
- else if(*copula==5) //Frank
- {
- /*double nu;
- t1 = -param[0]*u[j];
- t2 = -param[0]*v[j];
- t3 = exp(t1);
- t4 = exp(t2);
- t5 = 1-t3;
- t6 = 1-t4;
- nu = 1-exp(-param[0]);
- t7 = t5*t6;
- t8 = nu-t7;
- out[j] = -1/param[0]*log(t8/nu);*/
- if (param[0]>0) {
- t1=-log1p(exp(-param[0]) * expm1(param[0]-u[j]*param[0])/expm1(-param[0]));
- t2=-log1p(exp(-param[0]) * expm1(param[0]-v[j]*param[0])/expm1(-param[0]));
- out[j] = -log1mexp(t1+t2-log1mexp(param[0]))/param[0];
- } else {
- out[j] =-1/param[0] * log(1 + exp(-(-log((exp(-param[0] * u[j]) - 1)/(exp(-param[0]) - 1)) + -log((exp(-param[0] * v[j]) - 1)/(exp(-param[0]) - 1)))) * (exp(-param[0]) - 1));
- }
- }
- else if(*copula==6) //Joe
- {
- t1 = 1-u[j];
- t2 = 1-v[j];
- t3 = pow(t1,param[0]);
- t4 = pow(t2,param[0]);
- t5 = t3*t4;
- out[j] = 1-pow(t3+t4-t5,1/param[0]);
- }
- else if(*copula==7) //BB1
- {
- t1 = pow(u[j],-param[0]);
- t2 = pow(v[j],-param[0]);
- t3 = t1-1;
- t4 = t2-1;
- t5 = pow(t3,param[1]);
- t6 = pow(t4,param[1]);
- t7 = t5+t6;
- t8 = pow(t7,1/param[1]);
- out[j] = pow(1+t8,-1/param[0]);
- }
- else if(*copula==8) //BB6
- {
- t1 = 1-u[j];
- t2 = 1-v[j];
- t3 = pow(t1,param[0]);
- t4 = pow(t2,param[0]);
- t5 = 1-t3;
- t6 = 1-t4;
- t7 = -log(t5);
- t8 = -log(t6);
- t9 = pow(t7,param[1]);
- t10 = pow(t8,param[1]);
- t11 = t9+t10;
- t12 = pow(t11,1/param[1]);
- t13 = exp(-t12);
- t14 = 1-t13;
- out[j] = 1-pow(t14,1/param[0]);
- }
- else if(*copula==9) //BB7
- {
- t1 = 1-u[j];
- t2 = 1-v[j];
- t3 = pow(t1,param[0]);
- t4 = pow(t2,param[0]);
- t5 = 1-t3;
- t6 = 1-t4;
- t7 = pow(t5,-param[1]);
- t8 = pow(t6,-param[1]);
- t9 = t7+t8-1;
- t10 = pow(t9,-1/param[1]);
- t11 = 1-t10;
- t12 = pow(t11,1/param[0]);
- out[j] = 1-t12;
- }
- else if(*copula==10) //BB8
+ int j;
+ double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14;
+
+ for(j=0;j<*n;j++)
+ {
+ if(u[j]>UMAX && v[j]>UMAX){ out[j]=1;}
+ else if(u[j]>UMAX){ out[j]=v[j];}
+ else if(v[j]>UMAX){ out[j]=u[j];}
+ else if(u[j]<UMIN || v[j]<UMIN){ out[j]=0;}
+ else
{
- double nu;
- t1 = param[1]*u[j];
- t2 = param[1]*v[j];
- t3 = 1-t1;
- t4 = 1-t2;
- t5 = pow(t3,param[0]);
- t6 = pow(t4,param[0]);
- t7 = 1-t5;
- t8 = 1-t6;
- nu = 1-param[1];
- nu = pow(nu,param[0]);
- nu = 1-nu;
- nu = 1/nu;
- t9 = 1-nu*t7*t8;
- t10 = pow(t9,1/param[0]);
- out[j] = 1/param[1]*(1-t10);
+ if(*copula==3) //Clayton
+ {
+ t1 = pow(u[j],-param[0]);
+ t2 = pow(v[j],-param[0]);
+ t3 = t1+t2-1;
+ out[j] = pow(t3,-1/param[0]);
+ }
+ else if(*copula==4) //Gumbel
+ {
+ t1 = -log(u[j]);
+ t2 = -log(v[j]);
+ t3 = pow(t1,param[0]);
+ t4 = pow(t2,param[0]);
+ t5 = t3+t4;
+ t6 = -pow(t5,1/param[0]);
+ out[j] = exp(t6);
+ }
+ else if(*copula==5) //Frank
+ {
+ /*double nu;
+ t1 = -param[0]*u[j];
+ t2 = -param[0]*v[j];
+ t3 = exp(t1);
+ t4 = exp(t2);
+ t5 = 1-t3;
+ t6 = 1-t4;
+ nu = 1-exp(-param[0]);
+ t7 = t5*t6;
+ t8 = nu-t7;
+ out[j] = -1/param[0]*log(t8/nu);*/
+ if (param[0]>0) {
+ t1=-log1p(exp(-param[0]) * expm1(param[0]-u[j]*param[0])/expm1(-param[0]));
+ t2=-log1p(exp(-param[0]) * expm1(param[0]-v[j]*param[0])/expm1(-param[0]));
+ out[j] = -log1mexp(t1+t2-log1mexp(param[0]))/param[0];
+ } else {
+ out[j] =-1/param[0] * log(1 + exp(-(-log((exp(-param[0] * u[j]) - 1)/(exp(-param[0]) - 1)) + -log((exp(-param[0] * v[j]) - 1)/(exp(-param[0]) - 1)))) * (exp(-param[0]) - 1));
+ }
+ }
+ else if(*copula==6) //Joe
+ {
+ t1 = 1-u[j];
+ t2 = 1-v[j];
+ t3 = pow(t1,param[0]);
+ t4 = pow(t2,param[0]);
+ t5 = t3*t4;
+ out[j] = 1-pow(t3+t4-t5,1/param[0]);
+ }
+ else if(*copula==7) //BB1
+ {
+ t1 = pow(u[j],-param[0]);
+ t2 = pow(v[j],-param[0]);
+ t3 = t1-1;
+ t4 = t2-1;
+ t5 = pow(t3,param[1]);
+ t6 = pow(t4,param[1]);
+ t7 = t5+t6;
+ t8 = pow(t7,1/param[1]);
+ out[j] = pow(1+t8,-1/param[0]);
+ }
+ else if(*copula==8) //BB6
+ {
+ t1 = 1-u[j];
+ t2 = 1-v[j];
+ t3 = pow(t1,param[0]);
+ t4 = pow(t2,param[0]);
+ t5 = 1-t3;
+ t6 = 1-t4;
+ t7 = -log(t5);
+ t8 = -log(t6);
+ t9 = pow(t7,param[1]);
+ t10 = pow(t8,param[1]);
+ t11 = t9+t10;
+ t12 = pow(t11,1/param[1]);
+ t13 = exp(-t12);
+ t14 = 1-t13;
+ out[j] = 1-pow(t14,1/param[0]);
+ }
+ else if(*copula==9) //BB7
+ {
+ t1 = 1-u[j];
+ t2 = 1-v[j];
+ t3 = pow(t1,param[0]);
+ t4 = pow(t2,param[0]);
+ t5 = 1-t3;
+ t6 = 1-t4;
+ t7 = pow(t5,-param[1]);
+ t8 = pow(t6,-param[1]);
+ t9 = t7+t8-1;
+ t10 = pow(t9,-1/param[1]);
+ t11 = 1-t10;
+ t12 = pow(t11,1/param[0]);
+ out[j] = 1-t12;
+ }
+ else if(*copula==10) //BB8
+ {
+ double nu;
+ t1 = param[1]*u[j];
+ t2 = param[1]*v[j];
+ t3 = 1-t1;
+ t4 = 1-t2;
+ t5 = pow(t3,param[0]);
+ t6 = pow(t4,param[0]);
+ t7 = 1-t5;
+ t8 = 1-t6;
+ nu = 1-param[1];
+ nu = pow(nu,param[0]);
+ nu = 1-nu;
+ nu = 1/nu;
+ t9 = 1-nu*t7*t8;
+ t10 = pow(t9,1/param[0]);
+ out[j] = 1/param[1]*(1-t10);
+ }
+ else if(*copula==41)
+ {
+ t1=qgamma(1.0-u[j],param[0],1,1,0);
+ t2=qgamma(1.0-v[j],param[0],1,1,0);
+ t3=pow(pow(t1,param[0])+pow(t2,param[0]),(1.0/param[0]));
+ out[j]=1.0-pgamma(t3,param[0],1,1,0);
+ }
}
- else if(*copula==41)
- {
- t1=qgamma(1.0-u[j],param[0],1,1,0);
- t2=qgamma(1.0-v[j],param[0],1,1,0);
- t3=pow(pow(t1,param[0])+pow(t2,param[0]),(1.0/param[0]));
- out[j]=1.0-pgamma(t3,param[0],1,1,0);
- }
- }
- }
-
+ }
+
}
@@ -473,197 +473,197 @@
void dbb1(double* u, double* v, int* n, double* param, double* out)
{
- int i;
- double th, de;
- double t1, t2, t3, t16, t17, t38, t39, t4, t5, t6, t7, t9, t10, t12, t13, t20, t24, t25, t27, t29, t32, t33, t34, t36, t43, t59;
-
- th = param[0];
- de = param[1];
-
- for(i=0;i<*n;i++)
- {
- t1 = pow(u[i],(-th));
- t2 = t1-1.0;
- t3 = pow(t2,de);
- t16 = 1./u[i];
- t17 = 1./t2;
- t38 = t1*t16;
- t39 = t38*t17;
- t4 = pow(v[i],(-th));
- t5 = t4-1.0;
- t6 = pow(t5,de);
- t7 = t3+t6;
- t9 = pow(t7,(1./de));
- t10 = 1.0+t9;
- t12 = pow(t10,(-1./th));
- t13 = t12*t9;
- t20 = 1./t10;
- t24 = t9*t9;
- t25 = t12*t24;
- t27 = 1./v[i];
- t29 = 1./t5;
- t32 = t7*t7;
- t33 = 1./t32;
- t34 = t10*t10;
- t36 = t33/t34;
- t43 = t4*th;
- t59 = t43*t27*t29;
-
- out[i] = t25*t6*t27*t4*t29*t36*t3*t39-t13*t6*t43*t27*t29*t33*t3*t38*t17*t20+
- t13*t3*t38*t17*t33*t20*t6*de*t59+t25*t3*t39*t36*t6*t59;
- }
-
+ int i;
+ double th, de;
+ double t1, t2, t3, t16, t17, t38, t39, t4, t5, t6, t7, t9, t10, t12, t13, t20, t24, t25, t27, t29, t32, t33, t34, t36, t43, t59;
+
+ th = param[0];
+ de = param[1];
+
+ for(i=0;i<*n;i++)
+ {
+ t1 = pow(u[i],(-th));
+ t2 = t1-1.0;
+ t3 = pow(t2,de);
+ t16 = 1./u[i];
+ t17 = 1./t2;
+ t38 = t1*t16;
+ t39 = t38*t17;
+ t4 = pow(v[i],(-th));
+ t5 = t4-1.0;
+ t6 = pow(t5,de);
+ t7 = t3+t6;
+ t9 = pow(t7,(1./de));
+ t10 = 1.0+t9;
+ t12 = pow(t10,(-1./th));
+ t13 = t12*t9;
+ t20 = 1./t10;
+ t24 = t9*t9;
+ t25 = t12*t24;
+ t27 = 1./v[i];
+ t29 = 1./t5;
+ t32 = t7*t7;
+ t33 = 1./t32;
+ t34 = t10*t10;
+ t36 = t33/t34;
+ t43 = t4*th;
+ t59 = t43*t27*t29;
+
+ out[i] = t25*t6*t27*t4*t29*t36*t3*t39-t13*t6*t43*t27*t29*t33*t3*t38*t17*t20+
+ t13*t3*t38*t17*t33*t20*t6*de*t59+t25*t3*t39*t36*t6*t59;
+ }
+
}
void dbb6(double* u, double* v, int* n, double* param, double* out)
{
- int i;
- double th, de;
- double t1, t2, t3, t4, t5, t12, t16, t32, t38, t39, t40, t47, t50, t61, t90, t6, t7, t8, t9, t10, t11, t13, t14, t35, t36, t37, t42, t48, t53, t56, t57, t59, t78, t80, t87, t93;
-
- th = param[0];
- de = param[1];
-
- for(i=0;i<*n;i++)
- {
- t1 = 1.0-u[i];
- t2 = pow(t1,th);
- t3 = 1.0-t2;
- t4 = log(t3);
- t5 = pow(-t4,de);
- t12 = 1/de;
- t16 = 1/th;
- t32 = de-1.0;
- t38 = 2.0*de;
- t39 = -1.0+t38;
- t40 = pow(-t4,t39);
- t47 = 3.0*de-1.0;
- t50 = pow(-t4,t32);
- t61 = pow(-t4,t47);
- t90 = pow(-t4,t38);
- t6 = 1.0-v[i];
- t7 = pow(t6,th);
- t8 = 1.0-t7;
- t9 = log(t8);
- t10 = pow(-t9,de);
- t11 = t5+t10;
- t13 = pow(t11,t12);
- t14 = exp(-t13);
- //t15 = 1.0-t14;
- //t17 = pow(t15,t16);
- t35 = pow(t11,-2.0*t32*t12);
- t36 = t35*th;
- t37 = exp(t13);
- t42 = pow(-t9,t39);
- t48 = pow(-t9,t47);
- t53 = t13*de;
- t56 = pow(-t9,t32);
- t57 = t37*t50*t56;
- t59 = t13*th;
- t78 = t37-1.0;
- t80 = pow(t78*t14,t16);
- t87 = t78*t78;
- t93 = pow(-t9,t38);
-
- out[i] = (2.0*t36*t37*t40*t42+t36*t37*t48*t50+t53*th*t57-t59*t57+
- t36*t37*t61*t56-2.0*t35*t40*t42-t35*t61*t56-t53*th*t50*t56+t59*t50*t56-
- t35*t48*t50) *t80*t7*t2/t3/t8/t87/(t90+2.0*t5*t10+t93)/t1/t6;
- }
-
+ int i;
+ double th, de;
+ double t1, t2, t3, t4, t5, t12, t16, t32, t38, t39, t40, t47, t50, t61, t90, t6, t7, t8, t9, t10, t11, t13, t14, t35, t36, t37, t42, t48, t53, t56, t57, t59, t78, t80, t87, t93;
+
+ th = param[0];
+ de = param[1];
+
+ for(i=0;i<*n;i++)
+ {
+ t1 = 1.0-u[i];
+ t2 = pow(t1,th);
+ t3 = 1.0-t2;
+ t4 = log(t3);
+ t5 = pow(-t4,de);
+ t12 = 1/de;
+ t16 = 1/th;
+ t32 = de-1.0;
+ t38 = 2.0*de;
+ t39 = -1.0+t38;
+ t40 = pow(-t4,t39);
+ t47 = 3.0*de-1.0;
+ t50 = pow(-t4,t32);
+ t61 = pow(-t4,t47);
+ t90 = pow(-t4,t38);
+ t6 = 1.0-v[i];
+ t7 = pow(t6,th);
+ t8 = 1.0-t7;
+ t9 = log(t8);
+ t10 = pow(-t9,de);
+ t11 = t5+t10;
+ t13 = pow(t11,t12);
+ t14 = exp(-t13);
+ //t15 = 1.0-t14;
+ //t17 = pow(t15,t16);
+ t35 = pow(t11,-2.0*t32*t12);
+ t36 = t35*th;
+ t37 = exp(t13);
+ t42 = pow(-t9,t39);
+ t48 = pow(-t9,t47);
+ t53 = t13*de;
+ t56 = pow(-t9,t32);
+ t57 = t37*t50*t56;
+ t59 = t13*th;
+ t78 = t37-1.0;
+ t80 = pow(t78*t14,t16);
+ t87 = t78*t78;
+ t93 = pow(-t9,t38);
+
+ out[i] = (2.0*t36*t37*t40*t42+t36*t37*t48*t50+t53*th*t57-t59*t57+
+ t36*t37*t61*t56-2.0*t35*t40*t42-t35*t61*t56-t53*th*t50*t56+t59*t50*t56-
+ t35*t48*t50) *t80*t7*t2/t3/t8/t87/(t90+2.0*t5*t10+t93)/t1/t6;
+ }
+
}
void dbb7(double* u, double* v, int* n, double* param, double* out)
{
- int i;
- double th, de;
- double t1, t2, t3, t4, t5, t6, t7, t8, t9, t11, t12, t14, t15, t16, t18, t20, t23, t24, t25, t27, t30, t31, t32, t35, t37, t42, t54;
-
- th = param[0];
- de = param[1];
-
- for(i=0;i<*n;i++)
- {
- t1 = 1.0-u[i];
- t2 = pow(t1,th);
- t3 = 1.0-t2;
- t4 = pow(t3,-de);
- t5 = 1.0-v[i];
- t6 = pow(t5,th);
- t7 = 1.0-t6;
- t8 = pow(t7,-de);
- t9 = t4+t8-1.0;
- t11 = pow(t9,-1.0/de);
- t12 = 1.0-t11;
- t14 = pow(t12,1.0/th);
- t15 = t11*t11;
- t16 = t14*t15;
- t18 = 1./t5;
- t20 = 1./t7;
- t23 = t9*t9;
- t24 = 1./t23;
- t25 = t12*t12;
- t27 = t24/t25;
- t30 = t2/t1;
- t31 = 1./t3;
- t32 = t30*t31;
- t35 = t14*t11;
- t37 = t6*th;
- t42 = 1./t12;
- t54 = t37*t18*t20;
-
- out[i] = -t16*t8*t6*t18*t20*t27*t4*t32 + t35*t8*t37*t18*t20*t24*t4*t30*t31*t42+
- t35*t4*t30*t31*t24*t42*t8*de*t54+t16*t4*t32*t27*t8*t54;
- }
-
+ int i;
+ double th, de;
+ double t1, t2, t3, t4, t5, t6, t7, t8, t9, t11, t12, t14, t15, t16, t18, t20, t23, t24, t25, t27, t30, t31, t32, t35, t37, t42, t54;
+
+ th = param[0];
+ de = param[1];
+
+ for(i=0;i<*n;i++)
+ {
+ t1 = 1.0-u[i];
+ t2 = pow(t1,th);
+ t3 = 1.0-t2;
+ t4 = pow(t3,-de);
+ t5 = 1.0-v[i];
+ t6 = pow(t5,th);
+ t7 = 1.0-t6;
+ t8 = pow(t7,-de);
+ t9 = t4+t8-1.0;
+ t11 = pow(t9,-1.0/de);
+ t12 = 1.0-t11;
+ t14 = pow(t12,1.0/th);
+ t15 = t11*t11;
+ t16 = t14*t15;
+ t18 = 1./t5;
+ t20 = 1./t7;
+ t23 = t9*t9;
+ t24 = 1./t23;
+ t25 = t12*t12;
+ t27 = t24/t25;
+ t30 = t2/t1;
+ t31 = 1./t3;
+ t32 = t30*t31;
+ t35 = t14*t11;
+ t37 = t6*th;
+ t42 = 1./t12;
+ t54 = t37*t18*t20;
+
+ out[i] = -t16*t8*t6*t18*t20*t27*t4*t32 + t35*t8*t37*t18*t20*t24*t4*t30*t31*t42+
+ t35*t4*t30*t31*t24*t42*t8*de*t54+t16*t4*t32*t27*t8*t54;
+ }
+
}
void dbb8(double* u, double* v, int* n, double* param, double* out)
{
- int i;
- double th, de;
- double t2, t3, t12, t16, t6, t7, t10, t11, t33, t38, t39, t49, t59, t69, t25, t26, t29, t44, t45, t50, t54, t62, t67;
-
- th = param[0];
- de = param[1];
-
- for(i=0;i<*n;i++)
- {
- t2 = 1.0-de*u[i];
- t3 = pow(t2,th);
- //t4 = 1.0-t3;
- t10 = 1.0-de;
- t11 = pow(t10,th);
- t12 = 1.0-t11;
- //t13 = 1/t12;
- t16 = 1/th;
- t33 = th*t3;
- t38 = 2.0*th;
- t39 = pow(t10,t38);
- t49 = pow(t2,t38);
- t59 = pow(t10,3.0*th);
- t69 = t12*t12;
- t6 = 1.0-de*v[i];
- t7 = pow(t6,th);
- //t8 = 1.0-t7;
- //t15 = 1.0-(1.0-t3)*t8*t13;
- //t17 = pow(t15,t16);
- t25 = t3*t7;
- t26 = t11-t7-t3+t25;
- t29 = pow(-t26/t12,t16);
- t44 = pow(t6,t38);
- t45 = t3*t44;
- t50 = t49*t7;
- t54 = t49*t44;
- t62 = -2.0*t25*t11+t25-t33*t7+3.0*t33*t7*t11-3.0*t33*t7*t39+t25*t39+
- 2.0* t45*t11-t45*t39+2.0*t50*t11-t50*t39-2.0*t54*t11+t54*t39+t54-
- t50-t45+t33*t7*t59;
- t67 = t26*t26;
- out[i] = -de*t29*t62/t6/t2/t67/t69;
- }
-
+ int i;
+ double th, de;
+ double t2, t3, t12, t16, t6, t7, t10, t11, t33, t38, t39, t49, t59, t69, t25, t26, t29, t44, t45, t50, t54, t62, t67;
+
+ th = param[0];
+ de = param[1];
+
+ for(i=0;i<*n;i++)
+ {
+ t2 = 1.0-de*u[i];
+ t3 = pow(t2,th);
+ //t4 = 1.0-t3;
+ t10 = 1.0-de;
+ t11 = pow(t10,th);
+ t12 = 1.0-t11;
+ //t13 = 1/t12;
+ t16 = 1/th;
+ t33 = th*t3;
+ t38 = 2.0*th;
+ t39 = pow(t10,t38);
+ t49 = pow(t2,t38);
+ t59 = pow(t10,3.0*th);
+ t69 = t12*t12;
+ t6 = 1.0-de*v[i];
+ t7 = pow(t6,th);
+ //t8 = 1.0-t7;
+ //t15 = 1.0-(1.0-t3)*t8*t13;
+ //t17 = pow(t15,t16);
+ t25 = t3*t7;
+ t26 = t11-t7-t3+t25;
+ t29 = pow(-t26/t12,t16);
+ t44 = pow(t6,t38);
+ t45 = t3*t44;
+ t50 = t49*t7;
+ t54 = t49*t44;
+ t62 = -2.0*t25*t11+t25-t33*t7+3.0*t33*t7*t11-3.0*t33*t7*t39+t25*t39+
+ 2.0* t45*t11-t45*t39+2.0*t50*t11-t50*t39-2.0*t54*t11+t54*t39+t54-
+ t50-t45+t33*t7*t59;
+ t67 = t26*t26;
+ out[i] = -de*t29*t62/t6/t2/t67/t69;
+ }
+
}
@@ -681,153 +681,153 @@
//////////////////////////////////////////////////////////////
/*
-void LL_mod(int* family, int* n, double* u, double* v, double* theta, double* nu, double* loglik)
-{
- double* negv;
- double* negu;
- negv = (double *) malloc(*n*sizeof(double));
- negu = (double *) malloc(*n*sizeof(double));
- double ntheta, nnu;
- int nfamily;
- ntheta = -*theta;
- nnu = -*nu;
+ void LL_mod(int* family, int* n, double* u, double* v, double* theta, double* nu, double* loglik)
+ {
+ double* negv;
+ double* negu;
+ negv = (double *) malloc(*n*sizeof(double));
+ negu = (double *) malloc(*n*sizeof(double));
+ double ntheta, nnu;
+ int nfamily;
+ ntheta = -*theta;
+ nnu = -*nu;
+
+ for(int i=0;i<*n;i++)
+ {
+ if(u[i]<UMIN) u[i]=UMIN;
+ else if(u[i]>UMAX) u[i]=UMAX;
+ if(v[i]<UMIN) v[i]=UMIN;
+ else if(v[i]>UMAX) v[i]=UMAX;
+ }
+
+ if((*family==43))
+ {
+ nfamily=3;
+ if(*theta > 0){
+ ntheta=2*(*theta)/(1-*theta);
+ LL(&nfamily, n, u, v, &ntheta, nu, loglik);
+ }else{
+ ntheta=-2*(*theta)/(1+*theta);
+ for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
+ LL(&nfamily, n, u, negv, &ntheta, &nnu, loglik);
+ }
+ }else if((*family==44))
+ {
+ nfamily=4;
+ if(*theta > 0){
+ ntheta=1/(1-*theta);
+ LL(&nfamily, n, u, v, &ntheta, nu, loglik);
+ }else{
+ ntheta=1/(1+*theta);
+ for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
+ LL(&nfamily, n, u, negv, &ntheta, &nnu, loglik);
+ }
+ }else{
+ if((*family==23) | (*family==24) | (*family==26) | (*family==27) | (*family==28) | (*family==29) | (*family==30) | (*family==61)) // 90? rotated copulas
+ {
+ nfamily = (*family)-20;
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vinecopula -r 136
Mehr Informationen über die Mailingliste Vinecopula-commits