leto/math--gsl — Issue #38: Add Bspline example from Jurgen Rose
Author: leto
Created: 2011-10-10 15:14:11+00:00
State: open

``` perl
#!/usr/bin/perl -w

use strict;

use Math::GSL::BSpline qw/:all/;
use Math::GSL::Vector qw/:all/;
use Math::GSL::Matrix qw/:all/;
use Math::GSL::RNG qw/:all/;
use Math::GSL::Multifit qw /:all/;
use Math::GSL::Randist qw/:all/;
use Math::GSL::Statistics qw /:all/;

# number of data points to fit
my $N=200;

# number of fit coefficients
my $NCOEFFS=12;

# nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4
my $NBREAK=$NCOEFFS -2;

my $n = $N;
my $ncoeffs = $NCOEFFS;
my $nbreak = $NBREAK;

my ($chisq, $Rsq, $dof, $tss);

my $seed=5;
#my $r = Math::GSL::RNG->new;
my $r = Math::GSL::RNG->new($gsl_rng_default,$seed);

# allocate a cubic bspline workspace (k = 4)
my $bw = gsl_bspline_alloc(4, $nbreak);

my $B = Math::GSL::Vector->new($ncoeffs);
my $x = Math::GSL::Vector->new($n);
my $y = Math::GSL::Vector->new($n);
my $X = Math::GSL::Matrix->new($n, $ncoeffs);
my $c = Math::GSL::Vector->new($ncoeffs);
my $w = Math::GSL::Vector->new($n);
my $cov = Math::GSL::Matrix->new($ncoeffs, $ncoeffs);
my $mw = gsl_multifit_linear_alloc($n, $ncoeffs);

printf("#m=0,S=0\n");
# this is the data to be fitted
for (my $i = 0; $i < $n; ++$i) {
my $xi = (15.0 / ($N - 1)) * $i;
my $yi = cos($xi) * exp(-0.1 * $xi);

my $sigma = 0.1 * $yi;
my $dy = gsl_ran_gaussian($r->raw(), $sigma);
$yi += $dy;

$x->set([$i], [$xi]);
$y->set([$i], [$yi]);
$w->set([$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 (my $i = 0; $i < $n; ++$i) {
my ($xi) = $x->get([$i]);

# compute B_j(xi) for all j
gsl_bspline_eval($xi, $B->raw(), $bw);

# fill in row i of X
for (my $j = 0; $j < $ncoeffs; ++$j) {
my ($Bj) = $B->get([$j]);
gsl_matrix_set($X->raw(), $i, $j, $Bj);
}
}

# do the fit
$chisq=gsl_multifit_wlinear($X->raw(), $w->raw(), $y->raw(), $c->raw(), $cov->raw(), $mw);

$dof = $n - $ncoeffs;

my @w_data=$w->as_list();
my @y_data=$y->as_list();
my $y_size=$y->length();
$tss = gsl_stats_wtss(\@w_data, 1, \@y_data, 1, $y->length());
$Rsq = 1.0 - $chisq / $tss;

printf STDERR "chisq/dof = %e, Rsq = %f\n",$chisq / $dof, $Rsq;

# output the smoothed curve
{
my ($xi, $yi, $yerr);

printf("#m=1,S=0\n");
for (my $xi = 0.0; $xi < 15.0; $xi += 0.1) {
gsl_bspline_eval($xi, $B->raw(), $bw);
my $rc;
($rc,$yi,$yerr)=gsl_multifit_linear_est($B->raw(), $c->raw(), $cov->raw());
printf("%f %f\n", $xi, $yi);
}
}

$r->free();
gsl_bspline_free($bw);
gsl_vector_free($B->raw());
gsl_vector_free($x->raw());
gsl_vector_free($y->raw());
gsl_matrix_free($X->raw());
gsl_vector_free($c->raw());
gsl_vector_free($w->raw());
gsl_matrix_free($cov->raw());
gsl_multifit_linear_free($mw);

# The output can be plotted with gnu graph.
#
# $ ./a.out > bspline.dat
# chisq/dof = 1.118217e+00, Rsq = 0.989771
# $ graph -T ps -X x -Y y -x 0 15 -y -1 1.3 < bspline.dat > bspline.ps
```