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
108
109
110
111
112
113
114
115
116
|
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
/* number of data points to fit */
#define N 200
/* number of fit coefficients */
#define NCOEFFS 8
/* nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */
#define NBREAK (NCOEFFS - 2)
int
main (void)
{
const size_t n = N;
const size_t ncoeffs = NCOEFFS;
const size_t nbreak = NBREAK;
size_t i, j;
gsl_bspline_workspace *bw;
gsl_vector *B;
double dy;
gsl_rng *r;
gsl_vector *c, *w;
gsl_vector *x, *y;
gsl_matrix *X, *cov;
gsl_multifit_linear_workspace *mw;
double chisq;
gsl_rng_env_setup();
r = gsl_rng_alloc(gsl_rng_default);
/* allocate a cubic bspline workspace (k = 4) */
bw = gsl_bspline_alloc(4, nbreak);
B = gsl_vector_alloc(ncoeffs);
x = gsl_vector_alloc(n);
y = gsl_vector_alloc(n);
X = gsl_matrix_alloc(n, ncoeffs);
c = gsl_vector_alloc(ncoeffs);
w = gsl_vector_alloc(n);
cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
mw = gsl_multifit_linear_alloc(n, ncoeffs);
printf("#m=0,S=0\n");
/* this is the data to be fitted */
for (i = 0; i < n; ++i)
{
double sigma;
double xi = (15.0/(N-1)) * i;
double yi = cos(xi) * exp(-0.1 * xi);
sigma = 0.1;
dy = gsl_ran_gaussian(r, sigma);
yi += dy;
gsl_vector_set(x, i, xi);
gsl_vector_set(y, i, yi);
gsl_vector_set(w, i, 1.0 / (sigma*sigma));
printf("%f %f\n", xi, yi);
}
/* use uniform breakpoints on [0, 15] */
gsl_bspline_knots_uniform(0.0, 15.0, bw);
/* construct the fit matrix X */
for (i = 0; i < n; ++i)
{
double xi = gsl_vector_get(x, i);
/* compute B_j(xi) for all j */
gsl_bspline_eval(xi, B, bw);
/* fill in row i of X */
for (j = 0; j < ncoeffs; ++j)
{
double Bj = gsl_vector_get(B, j);
gsl_matrix_set(X, i, j, Bj);
}
}
/* do the fit */
gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw);
/* output the smoothed curve */
{
double xi, yi, yerr;
printf("#m=1,S=0\n");
for (xi = 0.0; xi < 15.0; xi += 0.1)
{
gsl_bspline_eval(xi, B, bw);
gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
printf("%f %f\n", xi, yi);
}
}
gsl_rng_free(r);
gsl_bspline_free(bw);
gsl_vector_free(B);
gsl_vector_free(x);
gsl_vector_free(y);
gsl_matrix_free(X);
gsl_vector_free(c);
gsl_vector_free(w);
gsl_matrix_free(cov);
gsl_multifit_linear_free(mw);
return 0;
} /* main() */
|
