[Vinecopula-commits] r108 - pkg/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Di Aug 4 14:42:29 CEST 2015
Author: tnagler
Date: 2015-08-04 14:42:28 +0200 (Tue, 04 Aug 2015)
New Revision: 108
Modified:
pkg/src/hfunc.c
Log:
- correct Hinv2 for Tawn families (must change type before call to Hinv)
Modified: pkg/src/hfunc.c
===================================================================
--- pkg/src/hfunc.c 2015-08-03 12:56:46 UTC (rev 107)
+++ pkg/src/hfunc.c 2015-08-04 12:42:28 UTC (rev 108)
@@ -1,14 +1,14 @@
/*
-** hfunc.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.
-**
-*/
+ ** hfunc.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/hfunc.h"
@@ -25,125 +25,125 @@
void pcondbb1(double* u, double* v, int* n, double* param, double* out)
{
- int i;
- double th, de;
- double t1, t2, t3, t16, t17, t4, t5, t6, t7, t9, t10, t12, t13, t20;
-
- th = param[0];
- de = param[1];
- for(i=0;i<*n;i++)
- {
- t1 = pow(u[i],-th);
- t2 = t1-1.;
- t3 = pow(t2,de);
- t16 = 1./u[i];
- t17 = 1./t2;
- t4 = pow(v[i],-th);
- t5 = t4-1.;
- 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;
- out[i] = t13*t3*t1*t16*t17/t7*t20;
- }
-
+ int i;
+ double th, de;
+ double t1, t2, t3, t16, t17, t4, t5, t6, t7, t9, t10, t12, t13, t20;
+
+ th = param[0];
+ de = param[1];
+ for(i=0;i<*n;i++)
+ {
+ t1 = pow(u[i],-th);
+ t2 = t1-1.;
+ t3 = pow(t2,de);
+ t16 = 1./u[i];
+ t17 = 1./t2;
+ t4 = pow(v[i],-th);
+ t5 = t4-1.;
+ 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;
+ out[i] = t13*t3*t1*t16*t17/t7*t20;
+ }
+
}
void pcondbb6(double* u, double* v, int* n, double* param, double* out)
{
- int i;
- double th, de;
- double t1, t2, t3, t4, t5, t12, t16, t6, t7, t8, t9, t10, t11, t13, t14, t15, t17;
-
- 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;
- 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);
-
- out[i] = -t17*t13*t5*t2/t1/t3/t4/t11*t14/t15;
- }
-
+ int i;
+ double th, de;
+ double t1, t2, t3, t4, t5, t12, t16, t6, t7, t8, t9, t10, t11, t13, t14, t15, t17;
+
+ 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;
+ 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);
+
+ out[i] = -t17*t13*t5*t2/t1/t3/t4/t11*t14/t15;
+ }
+
}
void pcondbb7(double* u, double* v, int* n, double* param, double* out)
{
- int i;
- double th, de;
- double t1, t2, t3, t4, t6, t8, t9, t11, t12, t14;
-
- th = param[0];
- de = param[1];
-
- for(i=0;i<*n;i++)
- {
- t1 = 1.0-u[i];
- t2 = pow(t1,1.0*th);
- t3 = 1.0-t2;
- t4 = pow(t3,-1.0*de);
- t6 = pow(1.0-v[i],1.0*th);
- t8 = pow(1.0-t6,-1.0*de);
- t9 = t4+t8-1.0;
- t11 = pow(t9,-1.0/de);
- t12 = 1.0-t11;
- t14 = pow(t12,1.0/th);
-
- out[i] = t14*t11*t4*t2/t1/t3/t9/t12;
- }
-
+ int i;
+ double th, de;
+ double t1, t2, t3, t4, t6, t8, t9, t11, t12, t14;
+
+ th = param[0];
+ de = param[1];
+
+ for(i=0;i<*n;i++)
+ {
+ t1 = 1.0-u[i];
+ t2 = pow(t1,1.0*th);
+ t3 = 1.0-t2;
+ t4 = pow(t3,-1.0*de);
+ t6 = pow(1.0-v[i],1.0*th);
+ t8 = pow(1.0-t6,-1.0*de);
+ t9 = t4+t8-1.0;
+ t11 = pow(t9,-1.0/de);
+ t12 = 1.0-t11;
+ t14 = pow(t12,1.0/th);
+
+ out[i] = t14*t11*t4*t2/t1/t3/t9/t12;
+ }
+
}
void pcondbb8(double* u, double* v, int* n, double* param, double* out)
{
- int i;
- double th, de;
- double t2, t3, t12, t16, t6, t7, t8, t10, t11, t13, t15, t17;
-
- th = param[0];
- de = param[1];
-
- for(i=0;i<*n;i++)
- {
- t2 = 1.0-de*u[i];
- t3 = pow(t2,th);
- t10 = 1.0-de;
- t11 = pow(t10,th);
- t12 = 1.0-t11;
- t13 = 1/t12;
- t16 = 1/th;
- 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);
-
- out[i] = t17*t3/t2*t8*t13/t15;
- }
-
+ int i;
+ double th, de;
+ double t2, t3, t12, t16, t6, t7, t8, t10, t11, t13, t15, t17;
+
+ th = param[0];
+ de = param[1];
+
+ for(i=0;i<*n;i++)
+ {
+ t2 = 1.0-de*u[i];
+ t3 = pow(t2,th);
+ t10 = 1.0-de;
+ t11 = pow(t10,th);
+ t12 = 1.0-t11;
+ t13 = 1/t12;
+ t16 = 1/th;
+ 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);
+
+ out[i] = t17*t3/t2*t8*t13/t15;
+ }
+
}
@@ -151,228 +151,228 @@
// i.e. Hfunc1 and Hfunc2
void Hfunc1(int* family,int* n,double* u,double* v,double* theta,double* nu,double* out)
{
- double *negv, *negu;
- negv = (double *) malloc(*n* sizeof(double));
- negu = (double *) malloc(*n*sizeof(double));
- double ntheta, nnu;
- int nfamily, j, T=1;
- 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);
- Hfunc(&nfamily, n, u, v, &ntheta, &nnu, out);
- }else{
- ntheta=-2*(*theta)/(1+*theta);
- for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
- Hfunc(&nfamily, n, u, negv, &ntheta, &nnu, out);
- }
- }else if((*family)==44)
- {
- nfamily=4;
- if(*theta > 0){
- ntheta=1/(1-*theta);
- Hfunc (&nfamily, n, u, v, &ntheta, &nnu, out);
- }else{
- ntheta=1/(1+*theta);
- for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
- Hfunc (&nfamily, n, u, negv, &ntheta, &nnu, out);
- }
- }else{
-
- if(((*family==23) | (*family==24) | (*family==26) | (*family==27) | (*family==28) | (*family==29) | (*family==30) | (*family==61) ))
+ double *negv, *negu;
+ negv = (double *) malloc(*n* sizeof(double));
+ negu = (double *) malloc(*n*sizeof(double));
+ double ntheta, nnu;
+ int nfamily, j, T=1;
+ ntheta = -*theta;
+ nnu = -*nu;
+
+ for(int i=0;i<*n;i++)
{
- nfamily=(*family)-20;
- for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
- Hfunc (&nfamily, n, u, negv, &ntheta, &nnu, out);
+ 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;
}
- else if(((*family==33) | (*family==34) | (*family==36) | (*family==37) | (*family==38) | (*family==39) | (*family==40) | (*family==71) ))
- {
- nfamily=(*family)-30;
- for (int i = 0; i < *n; ++i) {negu[i]=1 - u[i];}
- Hfunc(&nfamily, n, negu, v, &ntheta, &nnu, out);
- for (int i = 0; i < *n; i++) {out[i]=1-out[i];};
- }
- // u and v enter in wrong order from BiCopHfunc and have to be treated accordingly
- else if(*family==104)
- {
- double par3=1;
- dC_du(v,u,n,theta,nu,&par3,out);
- }
- else if(*family==114)
- {
- double par3=1;
- for(j=0;j<*n;j++)
- {
- negv[j]= 1-v[j];
- negu[j]= 1-u[j];
- dC_du(&negv[j],&negu[j],&T,theta,nu,&par3,&out[j]);
- out[j]= 1-out[j];
- }
- }
- else if(*family==124)
- {
- double par3=1;
- for(j=0;j<*n;j++)
- {
- negv[j]= 1-v[j];
- dC_du(&negv[j],&u[j],&T,&ntheta,nu,&par3,&out[j]);
- }
- }
- else if(*family==134)
- {
- double par3=1;
- for(j=0;j<*n;j++)
- {
- negu[j]= 1-u[j];
- dC_du(&v[j],&negu[j],&T,&ntheta,nu,&par3,&out[j]);
- out[j]=1-out[j];
-
- }
- }
- else if(*family==204)
- {
- double par3=*nu;
- double par2=1;
- dC_du(v,u,n,theta,&par2,&par3,out);
- }
- else if(*family==214)
- {
- double par3=*nu;
- double par2=1;
- for(j=0;j<*n;j++)
- {
- negv[j]= 1-v[j];
- negu[j]= 1-u[j];
- dC_du(&negv[j],&negu[j],&T,theta,&par2,&par3,&out[j]);
- out[j]= 1-out[j];
- }
- }
- else if(*family==224)
- {
- double par3=*nu;
- double par2=1;
- for(j=0;j<*n;j++)
- {
- negv[j]= 1-v[j];
- dC_du(&negv[j],&u[j],&T,&ntheta,&par2,&par3,&out[j]);
- }
- }
- else if(*family==234)
- {
- double par3=*nu;
- double par2=1;
- for(j=0;j<*n;j++)
- {
- negu[j]= 1-u[j];
- dC_du(&v[j],&negu[j],&T,&ntheta,&par2,&par3,&out[j]);
- out[j]=1-out[j];
- }
- }
- else {
- Hfunc (family, n, u, v, theta, nu, out);
- }
- }
- // ensure that results are in [0,1]
- for(int j=0; j <* n; ++j){out[j] = MIN(MAX(out[j], 0), 1);}
- free(negv);
- free(negu);
+
+ if((*family)==43)
+ {
+ nfamily=3;
+ if(*theta > 0){
+ ntheta=2*(*theta)/(1-*theta);
+ Hfunc(&nfamily, n, u, v, &ntheta, &nnu, out);
+ }else{
+ ntheta=-2*(*theta)/(1+*theta);
+ for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
+ Hfunc(&nfamily, n, u, negv, &ntheta, &nnu, out);
+ }
+ }else if((*family)==44)
+ {
+ nfamily=4;
+ if(*theta > 0){
+ ntheta=1/(1-*theta);
+ Hfunc (&nfamily, n, u, v, &ntheta, &nnu, out);
+ }else{
+ ntheta=1/(1+*theta);
+ for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
+ Hfunc (&nfamily, n, u, negv, &ntheta, &nnu, out);
+ }
+ }else{
+
+ if(((*family==23) | (*family==24) | (*family==26) | (*family==27) | (*family==28) | (*family==29) | (*family==30) | (*family==61) ))
+ {
+ nfamily=(*family)-20;
+ for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
+ Hfunc (&nfamily, n, u, negv, &ntheta, &nnu, out);
+ }
+ else if(((*family==33) | (*family==34) | (*family==36) | (*family==37) | (*family==38) | (*family==39) | (*family==40) | (*family==71) ))
+ {
+ nfamily=(*family)-30;
+ for (int i = 0; i < *n; ++i) {negu[i]=1 - u[i];}
+ Hfunc(&nfamily, n, negu, v, &ntheta, &nnu, out);
+ for (int i = 0; i < *n; i++) {out[i]=1-out[i];};
+ }
+ // u and v enter in wrong order from BiCopHfunc and have to be treated accordingly
+ else if(*family==104)
+ {
+ double par3=1;
+ dC_du(v,u,n,theta,nu,&par3,out);
+ }
+ else if(*family==114)
+ {
+ double par3=1;
+ for(j=0;j<*n;j++)
+ {
+ negv[j]= 1-v[j];
+ negu[j]= 1-u[j];
+ dC_du(&negv[j],&negu[j],&T,theta,nu,&par3,&out[j]);
+ out[j]= 1-out[j];
+ }
+ }
+ else if(*family==124)
+ {
+ double par3=1;
+ for(j=0;j<*n;j++)
+ {
+ negv[j]= 1-v[j];
+ dC_du(&negv[j],&u[j],&T,&ntheta,nu,&par3,&out[j]);
+ }
+ }
+ else if(*family==134)
+ {
+ double par3=1;
+ for(j=0;j<*n;j++)
+ {
+ negu[j]= 1-u[j];
+ dC_du(&v[j],&negu[j],&T,&ntheta,nu,&par3,&out[j]);
+ out[j]=1-out[j];
+
+ }
+ }
+ else if(*family==204)
+ {
+ double par3=*nu;
+ double par2=1;
+ dC_du(v,u,n,theta,&par2,&par3,out);
+ }
+ else if(*family==214)
+ {
+ double par3=*nu;
+ double par2=1;
+ for(j=0;j<*n;j++)
+ {
+ negv[j]= 1-v[j];
+ negu[j]= 1-u[j];
+ dC_du(&negv[j],&negu[j],&T,theta,&par2,&par3,&out[j]);
+ out[j]= 1-out[j];
+ }
+ }
+ else if(*family==224)
+ {
+ double par3=*nu;
+ double par2=1;
+ for(j=0;j<*n;j++)
+ {
+ negv[j]= 1-v[j];
+ dC_du(&negv[j],&u[j],&T,&ntheta,&par2,&par3,&out[j]);
+ }
+ }
+ else if(*family==234)
+ {
+ double par3=*nu;
+ double par2=1;
+ for(j=0;j<*n;j++)
+ {
+ negu[j]= 1-u[j];
+ dC_du(&v[j],&negu[j],&T,&ntheta,&par2,&par3,&out[j]);
+ out[j]=1-out[j];
+ }
+ }
+ else {
+ Hfunc (family, n, u, v, theta, nu, out);
+ }
+ }
+ // ensure that results are in [0,1]
+ for(int j=0; j <* n; ++j){out[j] = MIN(MAX(out[j], 0), 1);}
+ free(negv);
+ free(negu);
}
void Hfunc2(int* family,int* n,double* v,double* u,double* theta,double* nu,double* out)
{
- double *negv, *negu;
- negv = (double *) malloc(*n * sizeof(double));
- negu = (double *) malloc(*n * sizeof(double));
- double ntheta, nnu;
- int nfamily;
- ntheta = -*theta;
- nnu = -*nu;
-
+ double *negv, *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);
- Hfunc (&nfamily, n, v, u, &ntheta, &nnu, out);
- }else{
- ntheta=-2*(*theta)/(1+*theta);
- for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
- Hfunc(&nfamily, n, negv, u, &ntheta, &nnu, out);
- for (int i = 0; i < *n; i++) {out[i]=1-out[i];};
- }
- }else if((*family)==44)
- {
- nfamily=4;
- if(*theta > 0){
- ntheta=1/(1-*theta);
- Hfunc (&nfamily, n, v, u, &ntheta, &nnu, out);
- }else{
- ntheta=1/(1+*theta);
- for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
- Hfunc(&nfamily, n, negv, u, &ntheta, &nnu, out);
- for (int i = 0; i < *n; i++) {out[i]=1-out[i];}; }
- }else{
-
- if(((*family==23) | (*family==24) | (*family==26) | (*family==27) | (*family==28) | (*family==29) | (*family==30) | (*family==61) ))
{
- nfamily=(*family)-20;
- for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
- Hfunc(&nfamily, n, negv, u, &ntheta, &nnu, out);
- for (int i = 0; i < *n; i++) {out[i]=1-out[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;
}
- else if(((*family==33) | (*family==34) | (*family==36) | (*family==37) | (*family==38) | (*family==39) | (*family==40) | (*family==71) ))
- {
- nfamily=(*family)-30;
- for (int i = 0; i < *n; ++i) {negu[i]=1 - u[i];}
- Hfunc(&nfamily, n, v, negu, &ntheta, &nnu, out);
- }
- // else if(*family==104 | *family==204 | *family==114 | *family==214)
- // {
- // u und v vertauschen (Unsauber, aber so sollte es funktionieren in unserer bisherigen Notation)
- // Hfunc(family,n,u,v,theta,nu,out);
- // }
- else if((*family==124) | (*family==224))
- {
- nfamily=(*family)-20;
- for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
- Hfunc(&nfamily, n, negv, u, &ntheta, nu, out);
- for (int i = 0; i < *n; i++) {out[i]=1-out[i];};
- }
- else if((*family==134) | (*family==234))
- {
- nfamily=(*family)-30;
- for (int i = 0; i < *n; ++i) {negu[i]=1 - u[i];}
- Hfunc(&nfamily, n, v, negu, &ntheta, nu, out);
- }
- else
- {
- Hfunc(family, n, v, u, theta, nu, out);
+
+ if((*family)==43)
+ {
+ nfamily=3;
+ if(*theta > 0){
+ ntheta=2*(*theta)/(1-*theta);
+ Hfunc (&nfamily, n, v, u, &ntheta, &nnu, out);
+ }else{
+ ntheta=-2*(*theta)/(1+*theta);
+ for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
+ Hfunc(&nfamily, n, negv, u, &ntheta, &nnu, out);
+ for (int i = 0; i < *n; i++) {out[i]=1-out[i];};
+ }
+ }else if((*family)==44)
+ {
+ nfamily=4;
+ if(*theta > 0){
+ ntheta=1/(1-*theta);
+ Hfunc (&nfamily, n, v, u, &ntheta, &nnu, out);
+ }else{
+ ntheta=1/(1+*theta);
+ for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
+ Hfunc(&nfamily, n, negv, u, &ntheta, &nnu, out);
+ for (int i = 0; i < *n; i++) {out[i]=1-out[i];}; }
+ }else{
+
+ if(((*family==23) | (*family==24) | (*family==26) | (*family==27) | (*family==28) | (*family==29) | (*family==30) | (*family==61) ))
+ {
+ nfamily=(*family)-20;
+ for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
+ Hfunc(&nfamily, n, negv, u, &ntheta, &nnu, out);
+ for (int i = 0; i < *n; i++) {out[i]=1-out[i];};
+ }
+ else if(((*family==33) | (*family==34) | (*family==36) | (*family==37) | (*family==38) | (*family==39) | (*family==40) | (*family==71) ))
+ {
+ nfamily=(*family)-30;
+ for (int i = 0; i < *n; ++i) {negu[i]=1 - u[i];}
+ Hfunc(&nfamily, n, v, negu, &ntheta, &nnu, out);
+ }
+ // else if(*family==104 | *family==204 | *family==114 | *family==214)
+ // {
+ // u und v vertauschen (Unsauber, aber so sollte es funktionieren in unserer bisherigen Notation)
+ // Hfunc(family,n,u,v,theta,nu,out);
+ // }
+ else if((*family==124) | (*family==224))
+ {
+ nfamily=(*family)-20;
+ for (int i = 0; i < *n; ++i) {negv[i]=1 - v[i];}
+ Hfunc(&nfamily, n, negv, u, &ntheta, nu, out);
+ for (int i = 0; i < *n; i++) {out[i]=1-out[i];};
+ }
+ else if((*family==134) | (*family==234))
+ {
+ nfamily=(*family)-30;
+ for (int i = 0; i < *n; ++i) {negu[i]=1 - u[i];}
+ Hfunc(&nfamily, n, v, negu, &ntheta, nu, out);
+ }
+ else
+ {
+ Hfunc(family, n, v, u, theta, nu, out);
+ }
}
- }
- // ensure that results are in [0,1]
- for(int i=0; i < *n; ++i) {out[i] = MIN(MAX(out[i], 0), 1);}
- free(negv);
- free(negu);
+ // ensure that results are in [0,1]
+ for(int i=0; i < *n; ++i) {out[i] = MIN(MAX(out[i], 0), 1);}
+ free(negv);
+ free(negu);
}
@@ -390,534 +390,534 @@
//////////////////////////////////////////////////////////////
void Hfunc(int* family, int* n, double* u, double* v, double* theta, double* nu, double* out)
{
- int j;
- double *h;
- h = Calloc(*n,double);
- double x;
-
- /*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;
- }*/
-//Rprintf("family in Hfunc: %d\n", *family);
-//Rprintf("theta=par1 in Hfunc: %f\n", *theta);
-//Rprintf("nu=par2 in Hfunc: %f\n", *nu);
- for(j=0;j<*n;j++)
- {
- if((v[j]==0) | ( u[j]==0)) h[j] = 0;
- else if (v[j]==1) h[j] = u[j];
- else
- {
- if(*family==0) //independent
- {
- h[j] = u[j];
- }
- else if(*family==1) //gaussian
- {
- x = (qnorm(u[j],0.0,1.0,1,0) - *theta*qnorm(v[j],0.0,1.0,1,0))/sqrt(1.0-pow(*theta,2.0));
- if (isfinite(x))
- h[j] = pnorm(x,0.0,1.0,1,0);
- else if ((qnorm(u[j],0.0,1.0,1,0) - *theta*qnorm(v[j],0.0,1.0,1,0)) < 0)
- h[j] = 0;
- else
- h[j] = 1;
- }
- else if(*family==2) //student
- {
- double t1, t2, mu, sigma2;
- t1 = qt(u[j],*nu,1,0); t2 = qt(v[j],*nu,1,0); mu = *theta*t2; sigma2 = ((*nu+t2*t2)*(1.0-*theta*(*theta)))/(*nu+1.0);
- h[j] = pt((t1-mu)/sqrt(sigma2),*nu+1.0,1,0);
- }
- else if(*family==3) //clayton
- {
- if(*theta == 0) h[j] = u[j] ;
- if(*theta < XEPS) h[j] = u[j] ;
- else
- {
- x = pow(u[j],-*theta)+pow(v[j],-*theta)-1.0 ;
- h[j] = pow(v[j],-*theta-1.0)*pow(x,-1.0-1.0/(*theta));
- if(*theta < 0)
- {
- if(x < 0) h[j] = 0;
- }
- }
- }
- else if(*family==4) //gumbel
- {
- if(*theta == 1) h[j] = u[j] ;
- else
- {
- h[j] = -(exp(-pow(pow(-log(v[j]),*theta)+pow(-log(u[j]),*theta),1.0/(*theta)))*pow(pow(-log(v[j]),*theta)+pow(-log(u[j]),*theta),1.0/(*theta)-1.0)*pow(-log(v[j]),*theta))/(v[j]*log(v[j]));
- }
- }
- else if(*family==5) //frank
- {
- if(*theta==0) h[j]=u[j];
- else
- {
- h[j] = -(exp(*theta)*(exp(*theta*u[j])-1.0))/(exp(*theta*v[j]+*theta*u[j])-exp(*theta*v[j]+*theta)-exp(*theta*u[j]+*theta)+exp(*theta));
- }
- }
- else if(*family==6) //joe
- {
- if(*theta==1) h[j]=u[j];
- else
- {
- h[j] = pow(pow(1.0-u[j],*theta) + pow(1.0-v[j],*theta) - pow(1.0-u[j],*theta)*pow(1.0-v[j],*theta),1.0/(*theta)-1) * pow(1.0-v[j],*theta-1.0)*(1-pow(1-u[j],*theta));
- }
- }
- else if(*family==7) //BB1
- {
- double* param;
- param = Calloc(2,double);
- param[0]=*theta;
- param[1]=*nu;
- int T=1;
- if(*nu==1)
- {
- if(*theta==0) h[j]=u[j];
- else h[j]=pow(pow(u[j],-*theta)+pow(v[j],-*theta)-1,-1/(*theta)-1)*pow(v[j],-*theta-1);
- }
- else if(*theta==0)
- {
- h[j]=-(exp(-pow(pow(-log(v[j]),*nu)+pow(-log(u[j]),*nu),1.0/(*nu)))*pow(pow(-log(v[j]),*nu)+pow(-log(u[j]),*nu),1.0/(*nu)-1.0)*pow(-log(v[j]),*nu))/(v[j]*log(v[j]));
- }
- else
- {
- pcondbb1(&v[j],&u[j],&T,param,&h[j]);
- }
- Free(param);
- }
- else if(*family==8) //BB6
- {
- double* param;
- param = Calloc(2,double);
- param[0]=*theta;
- param[1]=*nu;
- int T=1;
- if(*theta==1)
- {
- if(*nu==1) h[j]=u[j];
- else h[j]=-(exp(-pow(pow(-log(v[j]),*nu)+pow(-log(u[j]),*nu),1.0/(*nu)))*pow(pow(-log(v[j]),*nu)+pow(-log(u[j]),*nu),1.0/(*nu)-1.0)*pow(-log(v[j]),*nu))/(v[j]*log(v[j]));
- }
- else if(*nu==1)
- {
- h[j]=pow(pow(1.0-u[j],*theta) + pow(1.0-v[j],*theta) - pow(1.0-u[j],*theta)*pow(1.0-v[j],*theta),1.0/(*theta)-1) * pow(1.0-v[j],*theta-1.0)*(1-pow(1-u[j],*theta));
- }
- else
- {
- pcondbb6(&v[j],&u[j],&T,param,&h[j]);
- }
- Free(param);
- }
- else if(*family==9) //BB7
- {
- double* param;
- param = Calloc(2,double);
- param[0]=*theta;
- param[1]=*nu;
- int T=1;
- if(*theta==1)
- {
- if(*nu==0) h[j]=u[j];
- else h[j]=pow(pow(u[j],-*nu)+pow(v[j],-*nu)-1,-1/(*nu)-1)*pow(v[j],-*nu-1);
- }
- else if(*nu==0)
- {
- h[j] = pow(pow(1.0-u[j],*theta) + pow(1.0-v[j],*theta) - pow(1.0-u[j],*theta)*pow(1.0-v[j],*theta),1.0/(*theta)-1) * pow(1.0-v[j],*theta-1.0)*(1-pow(1-u[j],*theta));
- }
- else
- {
- pcondbb7(&v[j],&u[j],&T,param,&h[j]);
- }
- Free(param);
- }
- else if(*family==10) //BB8
- {
- double* param;
- param = Calloc(2,double);
- param[0]=*theta;
- param[1]=*nu;
- int T=1;
- if(*nu==0)
- {
- h[j]=u[j];
- }
- else if(*nu==1)
- {
- if(*theta==1) h[j]=u[j];
- else h[j]=pow(pow(1.0-u[j],*theta) + pow(1.0-v[j],*theta) - pow(1.0-u[j],*theta)*pow(1.0-v[j],*theta),1.0/(*theta)-1) * pow(1.0-v[j],*theta-1.0)*(1-pow(1-u[j],*theta));
- }
- else
- {
- pcondbb8(&v[j],&u[j],&T,param,&h[j]);
- }
- Free(param);
- }
- else if(*family==13) //rotated clayton (180?)
- {
- if(*theta == 0) h[j] = u[j] ;
- if(*theta < XEPS) h[j] = u[j] ;
- else
- {
- u[j]=1-u[j];
- v[j]=1-v[j];
- x = pow(u[j],-*theta)+pow(v[j],-*theta)-1.0 ;
- h[j] = pow(v[j],-*theta-1.0)*pow(x,-1.0-1.0/(*theta)); // pow(v[j],-*theta-1.0)*pow(pow(u[j],-*theta)+pow(v[j],-*theta)-1.0,-1.0-1.0/(*theta));
- h[j]= 1-h[j];
- u[j]=1-u[j];
- v[j]=1-v[j];
- }
- }
- else if(*family==14) //rotated gumbel (180?)
- {
- v[j]= 1-v[j];
- u[j]= 1-u[j];
- h[j]= -(exp(-pow(pow(-log(v[j]),*theta)+pow(-log(u[j]),*theta),1.0/(*theta)))*pow(pow(-log(v[j]),*theta)+pow(-log(u[j]),*theta),1.0/(*theta)-1.0)*pow(-log(v[j]),*theta))/(v[j]* log(v[j]));
- h[j]= 1-h[j];
- u[j]=1-u[j];
- v[j]=1-v[j];
- }
- else if(*family==16)
- {
- v[j]= 1-v[j];
- u[j]= 1-u[j];
- h[j] = pow(pow(1.0-u[j],*theta) + pow(1.0-v[j],*theta) - pow(1.0-u[j],*theta)*pow(1.0-v[j],*theta),1.0/(*theta)-1) * pow(1.0-v[j],*theta-1.0)*(1-pow(1-u[j],*theta));
- h[j]= 1-h[j];
- u[j]=1-u[j];
- v[j]=1-v[j];
- }
- else if(*family==17) //rotated BB1
- {
- double* param;
- param = Calloc(2,double);
- param[0]=*theta;
- param[1]=*nu;
- int T=1;
- if(*nu==1)
- {
- if(*theta==0) h[j]=u[j];
- else
- {
- h[j]=pow(pow(1-u[j],-*theta)+pow(1-v[j],-*theta)-1,-1/(*theta)-1)*pow(1-v[j],-*theta-1);
- h[j]= 1-h[j];
- }
- }
- else if(*theta==0)
- {
- h[j]=-(exp(-pow(pow(-log(1-v[j]),*nu)+pow(-log(1-u[j]),*nu),1.0/(*nu)))*pow(pow(-log(1-v[j]),*nu)+pow(-log(1-u[j]),*nu),1.0/(*nu)-1.0)*pow(-log(1-v[j]),*nu))/((1-v[j])*log(1-v[j]));
- h[j]= 1-h[j];
- }
- else
- {
- v[j]= 1-v[j];
- u[j]= 1-u[j];
- pcondbb1(&v[j],&u[j],&T,param,&h[j]);
- u[j]=1-u[j];
- v[j]=1-v[j];
- h[j]= 1-h[j];
- }
- Free(param);
- }
- else if(*family==18) //rotated BB6
- {
- double* param;
- param = Calloc(2,double);
- param[0]=*theta;
- param[1]=*nu;
- int T=1;
- if(*theta==1)
- {
- if(*nu==1) h[j]=u[j];
- else
- {
- h[j]=-(exp(-pow(pow(-log(1-v[j]),*nu)+pow(-log(1-u[j]),*nu),1.0/(*nu)))*pow(pow(-log(1-v[j]),*nu)+pow(-log(1-u[j]),*nu),1.0/(*nu)-1.0)*pow(-log(1-v[j]),*nu))/((1-v[j])*log(1-v[j]));
- h[j]= 1-h[j];
- }
- }
- else if(*nu==1)
- {
- h[j]=pow(pow(u[j],*theta) + pow(v[j],*theta) - pow(u[j],*theta)*pow(v[j],*theta),1.0/(*theta)-1) * pow(v[j],*theta-1.0)*(1-pow(u[j],*theta));
- h[j]= 1-h[j];
- }
- else
- {
- v[j]= 1-v[j];
- u[j]= 1-u[j];
- pcondbb6(&v[j],&u[j],&T,param,&h[j]);
- u[j]=1-u[j];
- v[j]=1-v[j];
- h[j]= 1-h[j];
- }
- Free(param);
- }
- else if(*family==19) //rotated BB7
- {
- double* param;
- param = Calloc(2,double);
- param[0]=*theta;
- param[1]=*nu;
- int T=1;
- if(*theta==1)
- {
- if(*nu==0) h[j]=u[j];
- else{
- h[j]=pow(pow(1-u[j],-*nu)+pow(1-v[j],-*nu)-1,-1/(*nu)-1)*pow(1-v[j],-*nu-1);
- h[j]= 1-h[j];
- }
- }
- else if(*nu==0)
- {
- h[j] = pow(pow(u[j],*theta) + pow(v[j],*theta) - pow(u[j],*theta)*pow(v[j],*theta),1.0/(*theta)-1) * pow(v[j],*theta-1.0)*(1-pow(u[j],*theta));
- h[j]= 1-h[j];
- }
- else
- {
- v[j]= 1-v[j];
- u[j]= 1-u[j];
- pcondbb7(&v[j],&u[j],&T,param,&h[j]);
- u[j]=1-u[j];
- v[j]=1-v[j];
- h[j]= 1-h[j];
- }
- Free(param);
- }
- else if(*family==20) //rotated BB8
- {
- double* param;
- param = Calloc(2,double);
- param[0]=*theta;
- param[1]=*nu;
- int T=1;
- if(*nu==0)
- {
- h[j]=u[j];
- }
- else if(*nu==1)
- {
- if(*theta==1) h[j]=u[j];
- else{
- h[j]=pow(pow(u[j],*theta) + pow(v[j],*theta) - pow(u[j],*theta)*pow(v[j],*theta),1.0/(*theta)-1) * pow(v[j],*theta-1.0)*(1-pow(u[j],*theta));
- h[j]= 1-h[j];
- }
- }
- else
- {
- v[j]= 1-v[j];
- u[j]= 1-u[j];
- pcondbb8(&v[j],&u[j],&T,param,&h[j]);
- u[j]=1-u[j];
- v[j]=1-v[j];
- h[j]= 1-h[j];
- }
- Free(param);
- }
- else if(*family==41)
- {
- double t1,t2,t3;
- t1=qgamma(1.0-u[j],*theta,1,1,0);
- t2=qgamma(1.0-v[j],*theta,1,1,0);
- t3=pow(pow(t1,*theta)+pow(t2,*theta),(1.0/(*theta)));
- h[j]=exp(-t3+t1);
- }
- else if(*family==104)
- {
- int T=1;
- double par3=1;
- dC_dv(&u[j],&v[j],&T,theta,nu,&par3,&h[j]);
- }
- else if(*family==114)
- {
- int T=1;
- double par3=1;
- v[j]= 1-v[j];
- u[j]= 1-u[j];
- dC_dv(&u[j],&v[j],&T,theta,nu,&par3,&h[j]);
- u[j]=1-u[j];
- v[j]=1-v[j];
- h[j]= 1-h[j];
- }
- else if(*family==204)
- {
- int T=1;
- double par3=*nu, par2=1;
- dC_dv(&u[j],&v[j],&T,theta,&par2,&par3,&h[j]);
- }
- else if(*family==214)
- {
- int T=1;
- double par3=*nu, par2=1;
- v[j]= 1-v[j];
- u[j]= 1-u[j];
- dC_dv(&u[j],&v[j],&T,theta,&par2,&par3,&h[j]);
- u[j]=1-u[j];
- v[j]=1-v[j];
- h[j]= 1-h[j];
- }
+ int j;
+ double *h;
+ h = Calloc(*n,double);
+ double x;
+
+ /*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;
+ }*/
+ //Rprintf("family in Hfunc: %d\n", *family);
+ //Rprintf("theta=par1 in Hfunc: %f\n", *theta);
+ //Rprintf("nu=par2 in Hfunc: %f\n", *nu);
+ for(j=0;j<*n;j++)
+ {
+ if((v[j]==0) | ( u[j]==0)) h[j] = 0;
+ else if (v[j]==1) h[j] = u[j];
+ else
+ {
+ if(*family==0) //independent
+ {
+ h[j] = u[j];
+ }
+ else if(*family==1) //gaussian
+ {
+ x = (qnorm(u[j],0.0,1.0,1,0) - *theta*qnorm(v[j],0.0,1.0,1,0))/sqrt(1.0-pow(*theta,2.0));
+ if (isfinite(x))
+ h[j] = pnorm(x,0.0,1.0,1,0);
+ else if ((qnorm(u[j],0.0,1.0,1,0) - *theta*qnorm(v[j],0.0,1.0,1,0)) < 0)
+ h[j] = 0;
+ else
+ h[j] = 1;
+ }
+ else if(*family==2) //student
+ {
+ double t1, t2, mu, sigma2;
+ t1 = qt(u[j],*nu,1,0); t2 = qt(v[j],*nu,1,0); mu = *theta*t2; sigma2 = ((*nu+t2*t2)*(1.0-*theta*(*theta)))/(*nu+1.0);
+ h[j] = pt((t1-mu)/sqrt(sigma2),*nu+1.0,1,0);
+ }
+ else if(*family==3) //clayton
+ {
+ if(*theta == 0) h[j] = u[j] ;
+ if(*theta < XEPS) h[j] = u[j] ;
+ else
+ {
+ x = pow(u[j],-*theta)+pow(v[j],-*theta)-1.0 ;
+ h[j] = pow(v[j],-*theta-1.0)*pow(x,-1.0-1.0/(*theta));
+ if(*theta < 0)
+ {
+ if(x < 0) h[j] = 0;
+ }
+ }
+ }
+ else if(*family==4) //gumbel
+ {
+ if(*theta == 1) h[j] = u[j] ;
+ else
+ {
+ h[j] = -(exp(-pow(pow(-log(v[j]),*theta)+pow(-log(u[j]),*theta),1.0/(*theta)))*pow(pow(-log(v[j]),*theta)+pow(-log(u[j]),*theta),1.0/(*theta)-1.0)*pow(-log(v[j]),*theta))/(v[j]*log(v[j]));
+ }
+ }
+ else if(*family==5) //frank
+ {
+ if(*theta==0) h[j]=u[j];
+ else
+ {
+ h[j] = -(exp(*theta)*(exp(*theta*u[j])-1.0))/(exp(*theta*v[j]+*theta*u[j])-exp(*theta*v[j]+*theta)-exp(*theta*u[j]+*theta)+exp(*theta));
+ }
+ }
+ else if(*family==6) //joe
+ {
+ if(*theta==1) h[j]=u[j];
+ else
+ {
+ h[j] = pow(pow(1.0-u[j],*theta) + pow(1.0-v[j],*theta) - pow(1.0-u[j],*theta)*pow(1.0-v[j],*theta),1.0/(*theta)-1) * pow(1.0-v[j],*theta-1.0)*(1-pow(1-u[j],*theta));
+ }
+ }
+ else if(*family==7) //BB1
+ {
+ double* param;
+ param = Calloc(2,double);
+ param[0]=*theta;
+ param[1]=*nu;
+ int T=1;
+ if(*nu==1)
+ {
+ if(*theta==0) h[j]=u[j];
+ else h[j]=pow(pow(u[j],-*theta)+pow(v[j],-*theta)-1,-1/(*theta)-1)*pow(v[j],-*theta-1);
+ }
+ else if(*theta==0)
+ {
+ h[j]=-(exp(-pow(pow(-log(v[j]),*nu)+pow(-log(u[j]),*nu),1.0/(*nu)))*pow(pow(-log(v[j]),*nu)+pow(-log(u[j]),*nu),1.0/(*nu)-1.0)*pow(-log(v[j]),*nu))/(v[j]*log(v[j]));
+ }
+ else
+ {
+ pcondbb1(&v[j],&u[j],&T,param,&h[j]);
+ }
+ Free(param);
+ }
+ else if(*family==8) //BB6
+ {
+ double* param;
+ param = Calloc(2,double);
+ param[0]=*theta;
+ param[1]=*nu;
+ int T=1;
+ if(*theta==1)
+ {
+ if(*nu==1) h[j]=u[j];
+ else h[j]=-(exp(-pow(pow(-log(v[j]),*nu)+pow(-log(u[j]),*nu),1.0/(*nu)))*pow(pow(-log(v[j]),*nu)+pow(-log(u[j]),*nu),1.0/(*nu)-1.0)*pow(-log(v[j]),*nu))/(v[j]*log(v[j]));
+ }
+ else if(*nu==1)
+ {
+ h[j]=pow(pow(1.0-u[j],*theta) + pow(1.0-v[j],*theta) - pow(1.0-u[j],*theta)*pow(1.0-v[j],*theta),1.0/(*theta)-1) * pow(1.0-v[j],*theta-1.0)*(1-pow(1-u[j],*theta));
+ }
+ else
+ {
+ pcondbb6(&v[j],&u[j],&T,param,&h[j]);
+ }
+ Free(param);
+ }
+ else if(*family==9) //BB7
+ {
+ double* param;
+ param = Calloc(2,double);
+ param[0]=*theta;
+ param[1]=*nu;
+ int T=1;
+ if(*theta==1)
+ {
+ if(*nu==0) h[j]=u[j];
+ else h[j]=pow(pow(u[j],-*nu)+pow(v[j],-*nu)-1,-1/(*nu)-1)*pow(v[j],-*nu-1);
+ }
+ else if(*nu==0)
+ {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vinecopula -r 108
Mehr Informationen über die Mailingliste Vinecopula-commits