1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
| #include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sf_exp.h>
#include "cubature.h"
int count;
int n;
double bivar_cdf(unsigned dim, const double * x1, void *data)
{
int i,j;
int s;
double det,val, S;
double *R;
double P =dim;
R = (double *) malloc(1 * sizeof(double));
gsl_permutation* p = gsl_permutation_alloc(dim);
gsl_vector * x = gsl_vector_alloc (dim);
gsl_vector * xs = gsl_vector_alloc (dim);
gsl_matrix * A = gsl_matrix_alloc (dim, dim);
gsl_matrix * At = gsl_matrix_alloc (dim, dim);
gsl_matrix * SIG = gsl_matrix_alloc (dim, dim);
gsl_matrix *SIGi = gsl_matrix_alloc (dim, dim);
++count;
for (i=0; i < dim; i++){
for (j=0; j < dim; j++){
if (i < j) gsl_matrix_set (A,i,j,0);
else gsl_matrix_set (A,i,j,1);
}}
for (i=0; i < dim; i++) gsl_vector_set (x,i,x1[i]);
gsl_matrix_transpose_memcpy(At,A);
gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,1.0, A, At, 0.0, SIG);
gsl_linalg_LU_decomp( SIG, p, &s );
gsl_linalg_LU_invert(SIG,p,SIGi);
det = gsl_linalg_LU_det(SIG,s);
gsl_blas_dgemv(CblasNoTrans, 1.0, SIGi, x, 0.0, xs);
gsl_blas_ddot(x,xs,R);
S = exp ((-R[0]/2));
val =(1/(pow((2*M_PI),(P/2)))*(pow(abs(det),0.5)))*S;
return val;
}
int main()
{
int i;
double val, err;
double tol = 1.e-6;
double *xmin, *xmax;
unsigned dim = 6;
unsigned maxEval = 1000000;
printf("#==================================================================\n" );
printf("# Parameters:\n" );
printf("# dim = %d tol = %f maximum number of evals = %d\n",dim,tol,maxEval);
xmin = (double *) malloc(dim * sizeof(double));
xmax = (double *) malloc(dim * sizeof(double));
for (i = 0; i < dim; i++) {
xmin[i] = -99;
xmax[i] = 1.96;
}
count = 0;
adapt_integrate(bivar_cdf, 0, dim, xmin, xmax, maxEval, 0, tol, &val, &err);
printf("#------------------------------------------------------------------\n" );
printf("# Results:\n" );
printf("# Integration value = %g\n",val);
printf("# Est. err = %g\n", err);
printf("# Number of function evals required = %d\n", count);
printf("#==================================================================\n" );
free(xmax);
free(xmin);
return 0;
} |
