[Vinecopula-commits] r28 - in pkg/src: . include
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Do Okt 10 09:28:40 CEST 2013
Author: ulf
Date: 2013-10-10 09:28:40 +0200 (Thu, 10 Oct 2013)
New Revision: 28
Added:
pkg/src/evCopula.c
pkg/src/gof.c
pkg/src/include/evCopula.h
pkg/src/include/gof.h
pkg/src/include/pit.h
pkg/src/pit.c
Log:
Sorry, anscheinend hat das SVN die C-Dateien nicht commited
Added: pkg/src/evCopula.c
===================================================================
--- pkg/src/evCopula.c (rev 0)
+++ pkg/src/evCopula.c 2013-10-10 07:28:40 UTC (rev 28)
@@ -0,0 +1,265 @@
+#include "include/vine.h"
+#include "include/evCopula.h"
+#include <math.h>
+
+#define UMAX 1-1e-10
+
+#define UMIN 1e-10
+
+#define XEPS 1e-4
+
+// Some function for the Tawn copula
+
+// CDF
+void ta(double* t, int* n, double* par, double* par2, double* par3, double* out) //für CDF
+{
+ int i=0;
+ double t1,t2;
+ for(i=0; i<*n;i++)
+ {
+ t1=pow(*par2*t[i],*par);
+ t2=pow(*par3*(1.0-t[i]),*par);
+ out[i]=t1+t2;
+ }
+}
+
+//ta<-function(t,par,par2,par3) {(par2*t)^par+(par3*(1-t))^par}
+
+// Pickands A
+void Tawn(double* t, int* n, double* par, double* par2, double* par3, double* out) //für CDF
+{
+ int i=0, T=1;
+ double t1,t2,t3,t4;
+ for(i=0; i<*n;i++)
+ {
+ t1=(1.0-*par3)*(1.0-t[i]);
+ t2=(1.0-*par2)*t[i];
+ ta(t, &T, par, par2, par3, &t3);
+ t4=pow(t3,1.0/(*par));
+ out[i]=t1+t2+t4;
+ }
+}
+
+//Tawn<-function(t,par,par2,par3) {(1-par3)*(1-t)+(1-par2)*t+ta(t,par,par2,par3)^(1/par)}
+
+void TawnCDF(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out) // CDF-function
+{
+ int i=0, T=1;
+ double w, A;
+ for(i=0; i<*n;i++)
+ {
+ w=log(v[i])/log(u[i]*v[i]);
+ Tawn(&w, &T, par, par2, par3, &A); //!!!
+ out[i]=pow(u[i]*v[i],A);
+ }
+}
+
+
+//////////////////////////////////////////////////////////////////
+// PDF
+
+void ta2(double* t, int* n, double* par, double* par2, double* par3, double* out) //für PDF
+{
+ int i=0;
+ double t1,t2;
+ for(i=0; i<*n;i++)
+ {
+ t1=pow(*par3*t[i],*par);
+ t2=pow(*par2*(1.0-t[i]),*par);
+ out[i]=t1+t2;
+ }
+}
+
+void d1ta(double* t, int* n, double* par, double* par2, double* par3, double* out) //für PDF
+{
+ int i=0;
+ double t1,t2;
+ for(i=0; i<*n;i++)
+ {
+ t1=*par3 * pow((*par3*t[i]),*par-1.0);
+ t2=*par2 * pow(*par2*(1.0-t[i]),*par-1.0);
+ out[i]=*par*(t1-t2);
+ }
+}
+
+//d1ta<-function(t,par,par2,par3) {par*(par3*(par3*t)^(par-1)-par2*(par2*(1-t))^(par-1))}
+
+void d2ta(double* t, int* n, double* par, double* par2, double* par3, double* out) //für PDF
+{
+ int i=0;
+ double t1,t2;
+ for(i=0; i<*n;i++)
+ {
+ t1=pow(*par3,2) * pow(*par3*t[i],*par-2.0);
+ t2=pow(*par2,2) * pow(*par2*(1.0-t[i]),*par-2.0);
+ out[i]=*par*(*par-1) * (t1 + t2);
+ }
+}
+
+//d2ta<-function(t,par,par2,par3) {par*(par-1)*(par3^2*(par3*t)^(par-2)+par2^2*(par2*(1-t))^(par-2))}
+
+
+void Tawn2(double* t, int* n, double* par, double* par2, double* par3, double* out) //für PDF
+{
+ int i=0, T=1;
+ double t1,t2,t3,t4;
+ for(i=0; i<*n;i++)
+ {
+ t1=(1.0-*par2)*(1.0-t[i]);
+ t2=(1.0-*par3)*t[i];
+ ta2(&t[i], &T, par, par2, par3, &t3); //!!!
+ t4=pow(t3,1.0/(*par));
+ out[i]=t1+t2+t4;
+ }
+}
+
+//Tawn<-function(t,par,par2,par3) {(1-par2)*(1-t)+(1-par3)*t+ta(t,par,par2,par3)^(1/par)}
+
+void d1Tawn(double* t, int* n, double* par, double* par2, double* par3, double* out) //für PDF
+{
+ int i=0, T=1;
+ double t2,t1;
+ for(i=0; i<*n;i++)
+ {
+ ta2(&t[i], &T, par, par2, par3, &t1); //!!
+ d1ta(&t[i], &T, par, par2, par3, &t2);
+ out[i]=*par2-*par3+1.0/(*par) * pow(t1,(1.0/(*par)-1.0)) * t2;
+ }
+}
+
+//d1Tawn<-function(t,par,par2,par3) {par2-par3+1/par*ta(t,par,par2,par3)^(1/par-1)*d1ta(t,par,par2,par3)} Wie in Afunc2Deriv
+
+void d2Tawn(double* t, int* n, double* par, double* par2, double* par3, double* out) //für PDF
+{
+ int i=0, T=1;
+ double t2,t1,t3;
+ for(i=0; i<*n;i++)
+ {
+ ta2(&t[i], &T, par, par2, par3, &t1); //!!
+ d1ta(&t[i], &T, par, par2, par3, &t2);
+ d2ta(&t[i], &T, par, par2, par3, &t3);
+ out[i] = 1.0/(*par) * ( (1.0/(*par)-1.0) * pow(t1,(1.0/(*par)-2)) * pow(t2,2) + pow(t1,(1.0/(*par)-1)) * t3 );
+ }
+}
+
+//d2Tawn<-function(t,par,par2,par3) {1/par*((1/par-1)*ta(t,par,par2,par3)^(1/par-2)*d1ta(t,par,par2,par3)^2+ta(t,par,par2,par3)^(1/par-1)*d2ta(t,par,par2,par3))}
+
+// Ableitung von A nach u
+void dA_du(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out)
+{
+ int i=0, T=1;
+ double dA, dw, w;
+ for(i=0; i<*n;i++)
+ {
+ w=log(v[i]) / log(u[i]*v[i]);
+ dw=-log(v[i]) / (u[i]*pow(log(u[i]*v[i]),2.0));
+ d1Tawn(&w, &T, par, par2, par3, &dA);
+ out[i]=dA*dw;
+ }
+}
+
+//dA_du<-function(u,v,par,par2,par3) {evcBiCopAfuncDeriv(w(u,v),fam,par,par2,par3)*dw_du(u,v)}
+
+void dA_dv(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out)
+{
+ int i=0, T=1;
+ double dA, dw, w;
+ for(i=0; i<*n;i++)
+ {
+ w=log(v[i])/log(u[i]*v[i]);
+ dw=1.0 / (v[i]*log(u[i]*v[i])) - log(v[i]) / (v[i]*pow(log(u[i]*v[i]),2));
+ d1Tawn(&w, &T, par, par2, par3, &dA);
+ out[i]=dA*dw;
+ }
+}
+
+void dA_dudv(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out)
+{
+ int i=0, T=1;
+ double dA, dw_dv, dw_du, w, d2w_dudv, d2A;
+ for(i=0; i<*n;i++)
+ {
+ w=log(v[i])/log(u[i]*v[i]);
+ dw_du=-log(v[i])/(u[i]*pow(log(u[i]*v[i]),2));
+ dw_dv=1.0 / (v[i]*log(u[i]*v[i])) - log(v[i]) / (v[i]*pow(log(u[i]*v[i]),2));
+ d2w_dudv = 2*log(v[i]) / (v[i]*u[i]*pow(log(u[i]*v[i]),3)) - 1.0 / (v[i]*u[i]*pow(log(u[i]*v[i]),2));
+ d1Tawn(&w, &T, par, par2, par3, &dA);
+ d2Tawn(&w, &T, par, par2, par3, &d2A);
+ out[i]=d2A*dw_dv*dw_du + dA*d2w_dudv;
+ }
+}
+
+//d2A_dudv<-function(u,v,par,par2,par3) {evcBiCopAfunc2Deriv(w(u,v),fam,par,par2,par3)*dw_dv(u,v)*dw_du(u,v)+evcBiCopAfuncDeriv(w(u,v),fam,par,par2,par3)*d2w_dudv(u,v)}
+
+
+void TawnC(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out) // für PDF
+{
+ int i=0, T=1;
+ double w, A;
+ for(i=0; i<*n;i++)
+ {
+ w=log(v[i])/log(u[i]*v[i]);
+ Tawn2(&w, &T, par, par2, par3, &A); //!!!
+ out[i]=pow(u[i]*v[i],A);
+ }
+}
+
+// Ableitung von C nach u
+void dC_du(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out)
+{
+ int i=0, T=1;
+ double w, A, C, dA;
+ for(i=0; i<*n;i++)
+ {
+ w=log(v[i])/log(u[i]*v[i]);
+ Tawn2(&w, &T, par, par2, par3, &A); //!!
+ TawnC(&u[i], &v[i], &T, par, par2, par3, &C); //!!
+ dA_du(&u[i], &v[i], &T, par, par2, par3, &dA);
+ out[i]=C * (1.0 / u[i] * A + log(u[i]*v[i])*dA);
+ }
+}
+
+//dC_du<-function(u,v,par,par2,par3) {C(u,v,par,par2,par3) * (1/u*evcBiCopAfunc(w(u,v),fam,par,par2,par3)+log(u*v) * dA_du(u,v,par,par2,par3))}
+
+void TawnPDF(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out)
+{
+ int i=0, T=1;
+ double w, A, dC, t3, t4, t1, C, t5, t2;
+ for(i=0; i<*n;i++)
+ {
+ w=log(v[i])/log(u[i]*v[i]);
+ Tawn2(&w, n, par, par2, par3, &A);
+ dC_du(&u[i], &v[i], &T, par, par2, par3, &dC);
+ dA_du(&u[i], &v[i], &T, par, par2, par3, &t3);
+ dA_dv(&u[i], &v[i], &T, par, par2, par3, &t4);
+ t1=dC * (1.0/v[i] * A + log(u[i]*v[i]) * t4);
+ TawnC(&u[i], &v[i], &T, par, par2, par3, &C);
+ dA_dudv(&u[i], &v[i], &T, par, par2, par3, &t5);
+ t2=C * (1.0/u[i]*t4 + 1.0/v[i]*t3 + log(u[i]*v[i])*t5);
+ out[i]=t1+t2;
+ }
+}
+
+
+// d2C_dvdu<-function(u,v,par,par2,par3)
+// {
+// dC_du(u,v,par,par2,par3)*(1/v*evcBiCopAfunc(w(u,v),fam,par,par2,par3) + log(u*v)*dA_dv(u,v,par,par2,par3))+
+// C(u,v,par,par2,par3)*
+// (1/u*dA_dv(u,v,par,par2,par3) + 1/v*dA_du(u,v,par,par2,par3)+log(u*v)*d2A_dudv(u,v,par,par2,par3))
+// }
+
+
+// Ableitung von C nach v (fuer h-function)
+void dC_dv(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out)
+{
+ int i=0, T=1;
+ double w, A, C, dA;
+ for(i=0; i<*n;i++)
+ {
+ w=log(v[i])/log(u[i]*v[i]);
+ Tawn2(&w, &T, par, par2, par3, &A); //!!
+ TawnC(&u[i], &v[i], &T, par, par2, par3, &C); //!!
+ dA_dv(&u[i], &v[i], &T, par, par2, par3, &dA);
+ out[i]=C * (1.0 / v[i] * A + log(u[i]*v[i])*dA);
+ }
+}
\ No newline at end of file
Added: pkg/src/gof.c
===================================================================
--- pkg/src/gof.c (rev 0)
+++ pkg/src/gof.c 2013-10-10 07:28:40 UTC (rev 28)
@@ -0,0 +1,667 @@
+#include "include/vine.h"
+#include "include/memoryhandling.h"
+#include "include/gof.h"
+#include "include/rvinederiv2.h"
+#include "include/pit.h"
+#include "include/rvine.h"
+
+/////////////////////////////////////////////////////////////////////
+// Code form Daniel Berg, R-package copulaGOF
+// AD: Anderson-Darling GOF test
+// (Cumulative distribution function test)
+// INPUT:
+// cdf CDF for which to compute the test
+// n Length of cdf
+/////////////////////////////////////////////////////////////////////
+void ADtest(double* cdf, int* n, double* out)
+{
+ int j;
+ double sum=0.0;
+ for(j=0;j<*n;j++) sum += (2.0*((double)j+1.0)-1.0)*(log(cdf[j])+log(1.0-cdf[*n-1-j]));
+ *out = -(double)*n-(1.0/(double)*n)*sum;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+// Code form Daniel Berg, R-package copulaGOF
+// Function to compute cumulative distribution function of a uniform vector x ($\hat F(x)$)
+///////////////////////////////////////////////////////////////////////////////
+void CumDist(double* x, int* i_n, int* i_m, double* out)
+{
+ int i,j,n,m;
+ double *y;
+ n=*i_n; m=*i_m;
+ y = malloc(m*sizeof(double));
+ for(i=0;i<m;i++)
+ {
+ y[i]=0.0;
+ for(j=0;j<n;j++)
+ {
+ if(x[j]<=((double)i+1.0)/((double)m+1.0)) y[i] += 1.0/((double)n+1.0);
+ }
+ if(y[i]==0.0) y[i] = 1.0/((double)n+1.0);
+ out[i] = y[i];
+ }
+ free(y);
+}
+
+
+/////////////////////////////////////////////////////////////////////
+// Code form Daniel Berg, R-package copulaGOF
+// CvM: Cramer-von Mises GOF test
+// (Cumulative distribution function test)
+// INPUT:
+// cdf CDF for which to compute the test
+// n Length of cdf
+/////////////////////////////////////////////////////////////////////
+void CvMtest(double* cdf, int* n, double* out)
+{
+ int i;
+ double sum1=0.0,sum2=0.0;
+ for(i=0;i<*n;i++)
+ {
+ sum1 += pow(cdf[i],2.0);
+ sum2 += cdf[i]*(2.0*((double)i+1.0)+1.0);
+ }
+ *out = (double)*n/3.0 + (double)*n/((double)*n + 1.0)*sum1 - (double)*n/(pow((double)*n+1.0,2.0))*sum2;
+}
+
+
+/////////////////////////////////////////////////////////////////////
+// Code form Daniel Berg, R-package copulaGOF
+// KS: Kolmogorov-Smirnof GOF test
+// (Cumulative distribution function test)
+// INPUT:
+// cdf CDF for which to compute the test
+// n Length of cdf
+/////////////////////////////////////////////////////////////////////
+void KStest(double* cdf, int* n, double* out)
+{
+ int j;
+ double tmp, maxdist=0.0;
+ for(j=0;j<*n;j++)
+ {
+ tmp = MAX( fabs(cdf[j]-((double)j+1.0)/((double)*n+1.0)), fabs(cdf[j]-((double)j+2.0)/((double)*n+1.0)) );
+ if(tmp>maxdist) maxdist = tmp;
+ }
+ *out = sqrt((double)*n)*maxdist;
+}
+
+
+
+////////////////////////////////////////////////////////
+// Goodness-of-fit test based on White's information equality
+// by U. Schepsmeier
+///////////////////////////////////////////////////////////
+
+void White(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data, double* D, double* V)
+{
+ int i=0, dd=0, tt=0, k=1, j=0, kk=0, t=0, mm=0, dd2=0;
+ double *Dprime, *hess, *subhess, *der, *subder, *dat, *hess_red, *der_red;
+
+ for(i=0; i<(*d*(*d));i++)
+ {
+ if(family[i]!=0) dd++;
+ if(family[i]==2) tt++;
+ }
+ mm=(dd+tt)*(dd+tt+1)/2;
+ dd2=*d*(*d-1)/2;
+
+ //Allocate memory
+ //V = create_matrix((dd+tt)*(dd+tt+1)/2,(dd+tt)*(dd+tt+1)/2);
+ //D = malloc((dd+tt)*(dd+tt+1)/2*sizeof(double));
+ Dprime = malloc((dd+tt)*(dd+tt+1)/2*sizeof(double));
+ hess = malloc((dd2+tt)*(dd2+tt)*sizeof(double));
+ subhess = malloc((dd2+tt)*(dd2+tt)*sizeof(double));
+ der = malloc((dd2+tt)*(dd2+tt)*sizeof(double));
+ subder = malloc((dd2+tt)*(dd2+tt)*sizeof(double));
+ hess_red = malloc((dd+tt)*(dd+tt)*sizeof(double));
+ der_red = malloc((dd+tt)*(dd+tt)*sizeof(double));
+ dat = malloc(*d*sizeof(double));
+
+ // initialisieren
+ for(i=0;i<mm;i++)
+ {
+ Dprime[i]=0;
+ }
+ /*for(i=0;i<(dd+tt)*(dd+tt);i++)
+ {
+ hess[i]=0;
+ subhess[i]=0;
+ der[i]=0;
+ subder[i]=0;
+ }
+ for(i=0;i<(dd+tt)*(dd+tt+1)/2;i++)
+ {
+ for(j=0;j<(dd+tt)*(dd+tt+1)/2;j++)
+ {
+ V[i][j]=0;
+ }
+ }*/
+
+ for(t=0;t<*T;t++)
+ {
+ for(i=0; i<*d;i++)
+ {
+ dat[i]=data[(t+1)+(i*(*T))-1];
+ }
+ for(i=0;i<((dd2+tt)*(dd2+tt));i++)
+ {
+ hess[i]=0;
+ subhess[i]=0;
+ der[i]=0;
+ subder[i]=0;
+ }
+ hesse(&k, d, family, maxmat, matrix, condirect, conindirect, par, par2, dat, hess, subhess, der, subder);
+
+ // independence aus der Hesse herausnehmen
+ kk=0;
+ for(i=0;i<dd2+tt;i++)
+ {
+ for(j=0;j<dd2+tt;j++)
+ {
+ if(hess[i+1+((dd2+tt)*j)-1]!=0)
+ {
+ hess_red[kk]=hess[i+1+((dd2+tt)*j)-1];
+ der_red[kk]=der[i+1+((dd2+tt)*j)-1];
+ kk++;
+ }
+ }
+ }
+
+ kk=0;
+ for(i=0;i<(dd+tt);i++)
+ {
+ for(j=i;j<(dd+tt);j++)
+ {
+ Dprime[kk] = hess[(j+1)+((dd+tt)*i)-1] + der[(j+1)+((dd+tt)*i)-1];
+ D[kk] = D[kk] + (Dprime[kk]/(double)(*T));
+ kk++;
+ }
+ }
+
+ for(i=0;i<mm;i++)
+ {
+ for(j=0;j<mm;j++)
+ {
+ V[(i+1)+mm*j-1]+=(Dprime[i]*Dprime[j]/(double)(*T));
+ }
+ }
+ }
+
+ // Nicht fertig, da hier das Problem D%*%solve(V)%*%t(D) zu lösen ist
+
+ // Free memory
+ //free(D);
+ free(Dprime);
+ free(hess);
+ free(subhess);
+ free(der);
+ free(subder);
+ free(dat);
+ free(hess_red);
+ free(der_red);
+ //free_matrix(V,(dd+tt)*(dd+tt+1)/2);
+}
+
+
+///////////////////////////////////////////////////
+// Functions for PIT based GOF
+// by U. Schepsmeier
+///////////////////////////////////////////////////
+
+
+///////////////////////////////////////////////////////////////////////////////
+// Function to compute probability of observing rank i, given rank 1,...,i-1
+// Help function by Daniel Berg
+// Input: q (vector to be transformed, must be sorted ascending),
+// d(length of vector q), out(output)
+///////////////////////////////////////////////////////////////////////////////
+
+void ZStar(double* q, int* d, double* out)
+{
+ int i;
+ double *qprev;
+ qprev = malloc(*d*sizeof(double));
+ for(i=0;i<*d;i++)
+ {
+ if(i==0) qprev[i]=0.0;
+ else qprev[i] = q[i-1];
+ out[i] = 1.0 - pow((1.0-q[i])/(1.0-qprev[i]),*d-i);
+ if(out[i]==1.0) { out[i] = 0.9999999999; }
+ else if(out[i]==0.0) { out[i] = 1.0-0.9999999999; }
+ }
+ free(qprev);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Function to compare two numbers
+// Help function by Daniel Berg
+///////////////////////////////////////////////////////////////////////////////
+int comp_nums(double *num1, double *num2)
+{
+ if (*num1 < *num2) return -1;
+ else if (*num1 == *num2) return 0;
+ else return 1;
+}
+
+
+void Bj(int *T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* out, double* vv, double* vv2, int* calcupdate, int* method, int *alpha)
+{
+ int i=0, t=0, ii=0, j=0;
+ double *udata, **tmp, **u;
+ udata = malloc(*d*(*T)*sizeof(double));
+ tmp = create_matrix(*T,*d);
+ u = create_matrix(*T,*d);
+
+ RvinePIT(T, d, family, maxmat, matrix, condirect, conindirect, par, par2, data, udata, vv, vv2, calcupdate);
+
+ ii=0;
+ for(j=0;j<*T;j++)
+ {
+ if(*method==2 || *method==3)
+ {
+ for(i=0;i<*d;i++)
+ {
+ u[j][i] = udata[ii];
+ ii += 1;
+ }
+ qsort(u[j],*d,sizeof(double),(void *)comp_nums);
+ ZStar(u[j],d,tmp[j]); //Transformation von Berg and Bakken (2007)
+ }
+ else // Im Fall von Breymann ist es besser keine Transformation zu machen
+ {
+ for(i=0;i<*d;i++)
+ {
+ tmp[j][i] = udata[ii];
+ ii += 1;
+ }
+ }
+ }
+
+ for(t=0;t<*T;t++)
+ {
+ for(i=0;i<*d;i++)
+ {
+ if(*method==1)
+ tmp[t][i]=pow(qnorm(tmp[t][i],0.0,1.0,1,0),2);
+ else if(*method==2)
+ tmp[t][i]=fabs(tmp[t][i]-0.5);
+ else if(*method==3)
+ tmp[t][i]=pow(tmp[t][i]-0.5,*alpha);
+
+ out[t]+=tmp[t][i];
+ }
+ }
+
+ free(udata);
+ free_matrix(tmp,*T);
+ free_matrix(u,*T);
+}
+
+
+void SimulateBj(double* S, int *T, int* d, int* B, int* method, int *alpha, double* p)
+{
+ int i=0, t=0, m=0;
+ double *tmp, Sb=0, *ustar;
+ tmp = malloc(*d*sizeof(double));
+ ustar = malloc(*d*sizeof(double));
+
+ GetRNGstate();
+
+ for(t=0;t<*T;t++)
+ {
+ p[t]=0;
+ }
+
+ for(m=0;m<*B;m++)
+ {
+ for(i=1;i<=*d;i++) { tmp[i] = runif(0,1);}
+ qsort(tmp,*d,sizeof(double),(void *)comp_nums);
+ ZStar(tmp,d,ustar); //Transformation von Berg and Bakken (2007)
+
+ for(i=0;i<*d;i++)
+ {
+ if(*method==1)
+ tmp[i]=pow(qnorm(ustar[i],0.0,1.0,1,0),2);
+ else if(*method==2)
+ tmp[i]=fabs(ustar[i]-0.5);
+ else if(*method==3)
+ tmp[i]=pow(ustar[i]-0.5,*alpha);
+
+ Sb+=tmp[i];
+ }
+ for(t=0;t<*T;t++)
+ {
+ if(Sb<=S[t]) p[t]+=1.0/(1.0+(double)*B);
+ }
+ Sb=0;
+ }
+ for(t=0;t<*T;t++)
+ {
+ if(p[t] == 0) p[t] = 1.0/(1.0+(double)*B);
+ }
+
+ PutRNGstate();
+
+ free(tmp);
+ free(ustar);
+}
+
+
+void gofPIT_AD(int *T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* statistic, double* vv, double* vv2, int* calcupdate, int* method, int *alpha, int* B, int *statisticName)
+{
+ int t=0;
+ double *S, *helpvar, *Bhat;
+ S = malloc(*T*sizeof(double));
+ helpvar = malloc(*T*sizeof(double));
+ Bhat = malloc(*T*sizeof(double));
+
+ for(t=0;t<*T;t++)
+ {
+ S[t]=0;
+ helpvar[t]=0;
+ Bhat[t]=0;
+ }
+
+ Bj(T, d, family, maxmat, matrix, condirect, conindirect, par, par2, data, S, vv, vv2, calcupdate, method, alpha);
+
+ // Statistic berechnen
+ if(*B==0)
+ {
+ if(*method==1)
+ {
+ for(t=0;t<*T;t++)
+ {
+ Bhat[t]=pchisq(S[t],*d,1.0,0.0);
+ }
+ }
+ else
+ CumDist(S, T, T, Bhat);
+
+ if(*statisticName==1) //Anderson-Darling
+ ADtest(Bhat, T, statistic);
+ else if(*statisticName==2) //Kolmogorov-Smirnov
+ KStest(Bhat, T, statistic);
+ else if(*statisticName==3) //Cramer-von Mises
+ {
+ CvMtest(Bhat, T, statistic);
+ }
+
+ }
+ else //bootstrap
+ {
+ SimulateBj(S, T, d, B, method, alpha, helpvar);
+ CumDist(helpvar, T, T, Bhat);
+ if(*statisticName==1) //Anderson-Darling
+ ADtest(Bhat, T, statistic);
+ else if(*statisticName==2) //Kolmogorov-Smirnov
+ KStest(Bhat, T, statistic);
+ else if(*statisticName==3) //Cramer-von Mises
+ {
+ CvMtest(Bhat, T, statistic);
+ }
+ }
+
+ free(S);
+ free(helpvar);
+ free(Bhat);
+}
+
+
+void gofPIT_AD_pvalue(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* statistic, double* vv, double* vv2, int* calcupdate, int* method, int* alpha, int* B, double* pvalue, int *statisticName)
+{
+ int i=0, m=0, t=0, *f, B2=1000;
+ double *bdata, bstat=0;
+
+ f = malloc(*T*sizeof(int));
+ bdata = malloc(*d*(*T)*sizeof(double));
+
+ for(m=0;m<*B;m++)
+ {
+ MySample(T, T, f);
+ for(t=0;t<*T;t++)
+ {
+ for(i=0;i<*d;i++)
+ {
+ bdata[(t+1)+(*T*i)-1]=data[(f[t]+1)+(*T*i)-1];
+ }
+ }
+ bstat=0;
+ gofPIT_AD(T, d, family, maxmat, matrix, condirect, conindirect, par, par2, bdata,
+ &bstat, vv, vv2, calcupdate, method, alpha, &B2, statisticName);
+
+ if(bstat>=*statistic)
+ *pvalue+=1.0/(*B);
+ }
+
+ free(f);
+ free(bdata);
+}
+
+
+
+/* Equal probability sampling; with-replacement case */
+
+void MySample(int *k, int *n, int *y)
+{
+ int i;
+
+ GetRNGstate();
+ for (i = 0; i < *k; i++)
+ {
+ y[i] = (int) *n * unif_rand() + 1;
+ }
+ PutRNGstate();
+}
+
+
+////////////////////////////////////////////////////////////////
+
+// gof-test based on empirical copula process
+
+void gofECP(int* T, int* d, int* family, int* maxmat, int* matrix, int* conindirect, double* par, double* par2, double* data, double* statistic, int* statisticName)
+{
+ double *znull, *Chat1, *Chat2;
+ int T2=1000, i=0, t=0;
+ znull = malloc(*d*1000*sizeof(double));
+ Chat1 = malloc(*T*sizeof(double));
+ Chat2 = malloc(*T*sizeof(double));
+
+ for(t=0;t<T2;t++)
+ {
+ for(i=0;i<*d;i++)
+ {
+ znull[t+1+(T2*i)-1]=0;
+ }
+ }
+
+ SimulateRVine(&T2, d, family, maxmat, matrix, conindirect, par, par2, znull);
+
+
+ ChatZj(data, data, T, d, T, Chat1);
+ ChatZj(znull, data, T, d, &T2, Chat2);
+
+ *statistic=0;
+ if(*statisticName==3) //Cramer-von Mises test statistic
+ {
+ for(i=0;i<*T;i++)
+ {
+ *statistic+=pow(Chat1[i]-Chat2[i],2);
+ }
+ }
+ else if(*statisticName==2) // KS
+ {
+ for(i=0;i<*T;i++)
+ {
+ *statistic=MAX(fabs(Chat1[i]-Chat2[i]),*statistic);
+ }
+ *statistic=*statistic*sqrt(*T);
+ }
+
+ free(znull);
+ free(Chat1);
+ free(Chat2);
+}
+
+
+void gofECP_pvalue(int* T, int* d, int* family, int* maxmat, int* matrix, int* conindirect, double* par, double* par2, double* data, double* statistic, int* statisticName, double* pvalue, int* B)
+{
+ int i=0, m=0, t=0, *f;
+ double *bdata, bstat=0;
+
+ f = malloc(*T*sizeof(int));
+ bdata = malloc(*d*(*T)*sizeof(double));
+ //Rprintf("%f\n",*statistic);
+ for(m=0;m<*B;m++)
+ {
+ MySample(T, T, f);
+ for(t=0;t<*T;t++)
+ {
+ for(i=0;i<*d;i++)
+ {
+ bdata[(t+1)+(*T*i)-1]=data[(f[t]+1)+(*T*i)-1];
+ }
+ }
+ bstat=0;
+ gofECP(T, d, family, maxmat, matrix, conindirect, par, par2, bdata, &bstat, statisticName);
+ //Rprintf("%f ",bstat);
+ if(bstat>=*statistic)
+ *pvalue+=1.0/(*B);
+ }
+
+ free(f);
+ free(bdata);
+}
+
+
+// n = dim(u)[1]
+// m = dim(data)[1]
+// Chat vector of length n
+
+void ChatZj(double* data, double* u, int* n, int* d, int* m, double* Chat)
+{
+ int i,j,k;
+ double *helpvar;
+ helpvar=malloc(*m*sizeof(double));
+
+ for(j=0;j<*n;j++)
+ {
+ Chat[j]=0;
+ for(k=0;k<*m;k++)
+ {
+ helpvar[k]=0;
+ for(i=0;i<*d;i++)
+ {
+ if(data[k+1+(*m*i)-1]<=u[j+1+(*n*i)-1])
+ helpvar[k]++;
+ }
+ if(helpvar[k]==*d)
+ Chat[j]++;
+ }
+ Chat[j]/=(*m+1);
+ }
+
+ free(helpvar);
+}
+
+void C_ind(double* data, int* n, int* d, double* C)
+{
+ int t=0, i=0;
+
+ for(t=0;t<*n;t++)
+ {
+ for(i=0;i<*d;i++)
+ {
+ if(i==0)
+ C[t]=data[t+1+(*n*i)-1];
+ else
+ C[t]=C[t] * data[t+1+(*n*i)-1];
+ }
+
+ }
+}
+
+
+
+void gofECP2(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* vv, double* vv2, int* calcupdate, double* statistic, int* statisticName)
+{
+ double *udata, *Chat1, *Chat2;
+ int i=0, t=0;
+ udata = malloc(*d*(*T)*sizeof(double));
+ Chat1 = malloc(*T*sizeof(double));
+ Chat2 = malloc(*T*sizeof(double));
+
+ for(t=0;t<*T;t++)
+ {
+ for(i=0;i<*d;i++)
+ {
+ udata[t+1+(*T*i)-1]=0;
+ }
+ }
+ for(t=0;t<*T;t++)
+ {
+ Chat1[t]=0;
+ Chat2[t]=1;
+ }
+
+ RvinePIT(T, d, family, maxmat, matrix, condirect, conindirect, par, par2, data, udata, vv, vv2, calcupdate);
+ ChatZj(udata, udata, T, d, T, Chat1);
+
+ C_ind(udata,T,d,Chat2);
+
+ *statistic=0;
+ if(*statisticName==3) //Cramer-von Mises test statistic
+ {
+ for(i=0;i<*T;i++)
+ {
+ *statistic+=pow(Chat1[i]-Chat2[i],2);
+ }
+ }
+ else if(*statisticName==2) // KS
+ {
+ for(i=0;i<*T;i++)
+ {
+ *statistic=MAX(fabs(Chat1[i]-Chat2[i]),*statistic);
+ }
+ *statistic=*statistic*sqrt(*T);
+ }
+
+ free(udata);
+ free(Chat1);
+ free(Chat2);
+}
+
+void gofECP2_pvalue(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* vv, double* vv2, int* calcupdate, double* statistic, double* pvalue, int* statisticName, int* B)
+{
+ int i=0, m=0, t=0, *f;
+ double *bdata, bstat=0;
+
+ f = malloc(*T*sizeof(int));
+ bdata = malloc(*d*(*T)*sizeof(double));
+ //Rprintf("%f\n",*statistic);
+ for(m=0;m<*B;m++)
+ {
+ MySample(T, T, f);
+ for(t=0;t<*T;t++)
+ {
+ for(i=0;i<*d;i++)
+ {
+ bdata[(t+1)+(*T*i)-1]=data[(f[t]+1)+(*T*i)-1];
+ }
+ }
+ bstat=0;
+ gofECP2(T, d, family, maxmat, matrix, condirect, conindirect, par, par2, bdata, vv, vv2, calcupdate, &bstat, statisticName);
+ //Rprintf("%f ",bstat);
+ if(bstat>=*statistic)
+ *pvalue+=1.0/(*B);
+ }
+
+ free(f);
+ free(bdata);
+}
Added: pkg/src/include/evCopula.h
===================================================================
--- pkg/src/include/evCopula.h (rev 0)
+++ pkg/src/include/evCopula.h 2013-10-10 07:28:40 UTC (rev 28)
@@ -0,0 +1,31 @@
+#if !defined(EVCOPULA_H)
+#define EVCOPULA_H
+
+// Some function for the Tawn copula
+
+//CDF
+void ta(double* t, int* n, double* par, double* par2, double* par3, double* out);
+void Tawn(double* t, int* n, double* par, double* par2, double* par3, double* out);
+void TawnCDF(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out);
+
+///////////////////////////////////////////////////
+//PDF
+void ta2(double* t, int* n, double* par, double* par2, double* par3, double* out);
+void d1ta(double* t, int* n, double* par, double* par2, double* par3, double* out);
+void d2ta(double* t, int* n, double* par, double* par2, double* par3, double* out);
+
+void Tawn2(double* t, int* n, double* par, double* par2, double* par3, double* out);
+void d1Tawn(double* t, int* n, double* par, double* par2, double* par3, double* out);
+void d2Tawn(double* t, int* n, double* par, double* par2, double* par3, double* out);
+void dA_du(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out);
+void dA_dv(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out);
+void dA_dudv(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out);
+
+void TawnC(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out);
+void dC_du(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out);
+void TawnPDF(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out);
+
+// h-function
+void dC_dv(double* u, double* v, int* n, double* par, double* par2, double* par3, double* out);
+
+#endif
Added: pkg/src/include/gof.h
===================================================================
--- pkg/src/include/gof.h (rev 0)
+++ pkg/src/include/gof.h 2013-10-10 07:28:40 UTC (rev 28)
@@ -0,0 +1,31 @@
+#if !defined(GOF_H)
+#define GOF_H
+
+void ADtest(double* cdf, int* n, double* out); // Daniel Berg
+void CumDist(double* x, int* i_n, int* i_m, double* out); // Daniel Berg
+void CvMtest(double* cdf, int* n, double* out); // Daniel Berg
+void KStest(double* cdf, int* n, double* out); // Daniel Berg
+
+void White(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data, double* D, double* V);
+
+void Bj(int *T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* out, double* vv, double* vv2, int* calcupdate, int* method, int *alpha);
+void SimulateBj(double* S, int *T, int* d, int* B, int* method, int *alpha, double* p);
+void gofPIT_AD(int *T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* statistic, double* vv, double* vv2, int* calcupdate, int* method, int *alpha, int* B, int *statisticName);
+void gofPIT_AD_pvalue(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* statistic, double* vv, double* vv2, int* calcupdate, int* method, int* alpha, int* B, double* pvalue, int *statisticName);
+
+void MySample(int *k, int *n, int *y);
+
+void gofECP(int* T, int* d, int* family, int* maxmat, int* matrix, int* conindirect, double* par, double* par2, double* data, double* statistic, int* statisticName);
+void gofECP_pvalue(int* T, int* d, int* family, int* maxmat, int* matrix, int* conindirect, double* par, double* par2, double* data, double* statistic, int* statisticName, double* pvalue, int* B);
+void gofECP2(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* vv, double* vv2, int* calcupdate, double* statistic, int* statisticName);
+void gofECP2_pvalue(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* vv, double* vv2, int* calcupdate, double* statistic, double* pvalue,int* statisticName, int* B);
+
+void ChatZj(double* data, double* u, int* T, int* d, int* m, double* Chat);
+void C_ind(double* data, int* n, int* d, double* C);
+
+#endif
Added: pkg/src/include/pit.h
===================================================================
--- pkg/src/include/pit.h (rev 0)
+++ pkg/src/include/pit.h 2013-10-10 07:28:40 UTC (rev 28)
@@ -0,0 +1,5 @@
+// Probability integral transform
+void pit(int* n, int* d, int* family, int* type, double* par, double* nu, double* data, double* out);
+void RvinePIT(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* out, double* vv, double* vv2, int* calcupdate);
+
\ No newline at end of file
Added: pkg/src/pit.c
===================================================================
--- pkg/src/pit.c (rev 0)
+++ pkg/src/pit.c 2013-10-10 07:28:40 UTC (rev 28)
@@ -0,0 +1,262 @@
+//////////////////////////////////////////////////
+// PIT - Probability integral transform //
+// //
+// by Ulf Schepsmeier (2012) //
+//////////////////////////////////////////////////
+
+#include "include/vine.h"
+#include "include/memoryhandling.h"
+#include "include/likelihood.h"
+#include "include/pit.h"
+#include "include/hfunc.h"
+
+#define UMAX 1-1e-10
+
+#define UMIN 1e-10
+
+#define XEPS 1e-4
+
+//////////////////////////////////////////////////////////////
+// Function to transform a pair-copula construction (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)
+// par parameter values (at least d*(d-1)/2 parameters)
+////////////////////////////////////////////////////////////////
+
+void pit(int* T, int* d, int* family, int* type, double* par, double* nu, double* data, double* out)
+{
+ int i, j, in=1, k, **fam, tt;
+ double **v, **theta, **z, **ny, **x;
+
+ x = create_matrix(*d+1,*T);
+ v = create_matrix(*d+1,2*(*d)-1);
+ theta = create_matrix(*d+1,*d+1);
+ z = create_matrix(*d+1,*T);
+ ny = create_matrix(*d+1,*d+1);
+ fam = create_intmatrix(*d+1,*d+1);
+
+ k = 0;
+ for(i=0;i<*d;i++)
+ {
+ for (tt=0;tt<=*T-1;tt++ )
+ {
+ x[i][tt] = data[k];
+ k++;
+ }
+ }
+
+ //Initialize dependency parameters
+ k = 0;
+ for(i=0;i<*d-1;i++)
+ {
+ for(j=0;j<(*d-i-1);j++)
+ {
+ fam[i][j] = family[k];
+ ny[i][j] = nu[k];
+ theta[i][j] = par[k];
+ k++;
+ }
+ }
+
+ // Transform
+ if(*type==1) //Canonical vine
+ {
+ for(j=0;j<*T;j++)
+ {
+ z[0][j] = x[0][j];
+ for(i=1;i<*d;i++)
+ {
+ z[i][j]=x[i][j];
+ for(k=0;k<=(i-1);k++)
+ {
+ Hfunc1(&fam[k][i-k-1],&in, &z[i][j],&z[k][j],&theta[k][i-k-1],&ny[k][i-k-1],&z[i][j]);
+ }
+ }
+ }
+ }
+ else if(*type==2) //D-vine
+ {
+ for(j=0;j<*T;j++)
+ {
+ z[0][j] = x[0][j];
+ Hfunc1(&fam[0][0],&in, &x[1][j],&x[0][j],&theta[0][0],&ny[0][0],&z[1][j]);
+ v[1][0] = x[1][j];
+ Hfunc2(&fam[0][0],&in, &x[0][j],&x[1][j],&theta[0][0],&ny[0][0],&v[1][1]);
+ for(i=2;i<*d;i++)
+ {
+ Hfunc1(&fam[0][i-1],&in, &x[i][j],&x[i-1][j],&theta[0][i-1],&ny[0][i-1],&z[i][j]);
+ for(k=1;k<=(i-1);k++)
+ {
+ Hfunc1(&fam[k][i-k-1],&in, &z[i][j],&v[i-1][2*(k-1)+1],&theta[k][i-k-1],&ny[k][i-k-1],&z[i][j]);
+ }
+ if(i==(*d-1))
+ break;
+
+ v[i][0] = x[i][j];
+ Hfunc2(&fam[0][i-1],&in, &v[i-1][0],&v[i][0],&theta[0][i-1],&ny[0][i-1],&v[i][1]);
+ Hfunc1(&fam[0][i-1],&in, &v[i][0],&v[i-1][0],&theta[0][i-1],&ny[0][i-1],&v[i][2]);
+ if(i>2)
+ {
+ for(k=0;k<=(i-3);k++)
+ {
+ Hfunc2(&fam[k+1][i-k-2],&in, &v[i-1][2*k+1],&v[i][2*k+2],&theta[k+1][i-k-2],&ny[k+1][i-k-2],&v[i][2*k+3]);
+ Hfunc1(&fam[k+1][i-k-2],&in, &v[i][2*k+2],&v[i-1][2*k+1],&theta[k+1][i-k-2],&ny[k+1][i-k-2],&v[i][2*k+4]);
+ }
+ }
+ Hfunc2(&fam[i-1][0],&in, &v[i-1][2*i-3],&v[i][2*i-2],&theta[i-1][0],&ny[i-1][0],&v[i][2*i-1]);
+ }
+ }
+ }
+
+ //Write to output vector:
+ k = 0;
+ for(i=0;i<*d;i++)
+ {
+ for(j=0;j<*T;j++)
+ {
+ out[k] = z[i][j];
+ k ++;
+ }
+ }
+
+ //Free memory:
+ free_matrix(x,*d+1); free_matrix(v,*d+1); free_matrix(theta,*d+1); free_matrix(ny,*d+1); free_intmatrix(fam,*d+1); free_matrix(z,*d+1);
+}
+
+
+
+void RvinePIT(int* T, int* d, int* family, int* maxmat, int* matrix, int* condirect, int* conindirect, double* par, double* par2, double* data,
+ double* out, double* vv, double* vv2, int* calcupdate)
+{
+ int i, j, k, t, m, **fam;
+ double **x, **theta, **nu, ***vdirect, ***vindirect, **z;
+
+ //Allocate memory
+ x = create_matrix(*d,*T);
+ vdirect = create_3darray(*d,*d,*T);
+ vindirect = create_3darray(*d,*d,*T);
+ theta=create_matrix(*d,*d);
+ nu=create_matrix(*d,*d);
+ fam=create_intmatrix(*d,*d);
+ z = create_matrix(*d,*T);
+
+ //Initialize
+ k=0;
+ for(i=0;i<(*d);i++)
+ {
+ for (t=0;t<*T;t++ )
+ {
+ x[i][t] = data[k];
+ k++;
+ }
+ }
+
+ k=0;
+ for(i=0;i<(*d);i++)
+ {
+ for(j=0;j<(*d);j++)
+ {
+ theta[i][j]=par[(i+1)+(*d)*j-1] ;
+ nu[i][j]=par2[(i+1)+(*d)*j-1] ;
+ fam[i][j]=family[(i+1)+(*d)*j-1] ;
+ 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];
+ }
+ }
+ }
+
+ for(i=0;i<(*d);i++)
+ {
+ for(t=0;t<*T;t++ )
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vinecopula -r 28
Mehr Informationen über die Mailingliste Vinecopula-commits