Members/chinone/覚書

GSL

コンパイル

gcc -o hoge hoge.c -lgsl -lgslcblas -lm

乱数分布

/Members/chinone/覚書/Statistics

Polynomial Fit

データ作成

filemake_data.c
Everything is expanded.Everything is shortened.
  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
 
 
 
 
 
 
 
 
 
 
 
-
|
|
|
|
|
|
|
|
|
|
|
-
|
|
!
|
|
!
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_randist.h>
 
#define C0        (1.0)
#define C1        (3.0)
#define C2        (0.5)
 
#define SIGMA     (0.3)
 
int main(void)
{
  double x, y ;
  const gsl_rng_type *T ;
  gsl_rng *r ;
 
  gsl_rng_env_setup() ;
 
  T = gsl_rng_default ;
  r = gsl_rng_alloc(T) ;
 
 
 
  for( x=-2.0; x<2.1; x+=0.1){
    y = C0 + C1*cos(x) + C2*pow( x, 3.0) + gsl_ran_gaussian( r, SIGMA) ;
    printf("%10.5lf   %10.5lf\n", x, y) ;
  }
 
  return 0 ;
}
./make_data > test_data.dat

Fitting

filepolynomial_fit.c
Everything is expanded.Everything is shortened.
  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
117
118
119
120
121
122
123
124
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
-
|
|
|
|
|
|
|
|
|
|
|
-
|
|
!
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
-
|
|
|
|
!
|
|
|
|
|
-
|
|
|
|
|
|
!
|
|
|
|
|
|
|
-
|
|
|
!
|
|
|
|
|
|
|
|
|
-
|
|
|
|
!
|
|
|
|
!
 
 
 
-
|
|
-
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
!
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_multifit.h>
 
#define BUFFER_SIZE      (256)
#define MAX              (256)
#define PARAMETER_NUMBER (3)
#define SIGMA            (0.3)
 
#define C(i)             (gsl_vector_get(c,(i)))
#define COV(i,j)         (gsl_matrix_get(cov,(i),(j)))
 
double f(const int,const double) ;
 
int main(void)
{
  //*----------------- Input Data --------------------//
  FILE *fptr ;
  char buffer[BUFFER_SIZE] ;
  static double x[MAX], y[MAX] ;
  int i, j ;
  static int data_number ;
 
  if( (fptr=fopen( "./test_data.dat", "r")) == NULL )
    exit(EXIT_FAILURE) ;
 
  i = 0 ;
  while( fgets( buffer, BUFFER_SIZE, fptr)!=NULL ){
    sscanf( buffer, "%lf %lf", &x[i], &y[i]) ;
    i++ ;
  }
  data_number = i ;
  fclose(fptr) ;
 
  for( i=0; i<data_number; i++)
    printf("%10.5lf   %10.5lf\n", x[i], y[i]) ;
  printf("\n\n") ;
  //*-----------------------------------------------//
 
 
  //*------------------- vector & matrix ------------------------*//
  gsl_matrix *X, *cov ;
  gsl_vector *Y, *W, *c ;
 
  X = gsl_matrix_alloc( data_number, PARAMETER_NUMBER) ;
  Y = gsl_vector_alloc( data_number ) ;
  W = gsl_vector_alloc( data_number ) ;
  c = gsl_vector_alloc( PARAMETER_NUMBER ) ;
  cov = gsl_matrix_alloc( PARAMETER_NUMBER, PARAMETER_NUMBER) ;
 
  for( i=0; i<data_number; i++){
    gsl_vector_set( Y, i, y[i]) ;
    gsl_vector_set( W, i, 1.0/pow( SIGMA, 2.0)) ;
    for( j=0; j<PARAMETER_NUMBER; j++ )
      gsl_matrix_set( X, i, j, f(j, x[i])) ;
  }
  //*-----------------------------------------------------------*//
 
 
  //*------------------------- Polynomial Fit --------------------------*//
  double chisq ;
  {
    gsl_multifit_linear_workspace *work ;
    work = gsl_multifit_linear_alloc( data_number, PARAMETER_NUMBER) ;
 
    gsl_multifit_wlinear( X, W, Y, c, cov, &chisq, work) ;
 
    gsl_multifit_linear_free(work) ;
  }
  //*------------------------------------------------------------------*//
 
 
  //*----------------------- Print results ----------------------------*//
  for( j=0; j<PARAMETER_NUMBER; j++)
    printf("# c%d=%10.5lf\n", j, C(j)) ;
 
  for( i=0; i<PARAMETER_NUMBER; i++){
    for( j=0; j<PARAMETER_NUMBER; j++)
      printf("# %10.5lf  ", COV(i,j)) ;
    printf("\n") ;
  }
 
  printf("# chisq=%10.5lf\n", chisq) ;
  printf("# DOF: data_number-parameter_number=%d-%d=%d\n", data_number, PARAMETER_NUMBER, data_number-PARAMETER_NUMBER) ;
  printf("# reduced chisq=%10.5lf\n", chisq/(data_number-PARAMETER_NUMBER)) ;
  //*------------------------------------------------------------------*//
 
 
  //*----------------------- Output Fitted function ------------------*//
  double tmp ;
  for( i=0; i<data_number; i++){
    tmp = 0.0 ;
    for( j=0; j<PARAMETER_NUMBER; j++)
      tmp += C(j)*f(j, x[i]) ;
    printf("%10.5lf   %10.5lf\n", x[i], tmp) ;
  }
  //*-----------------------------------------------------------------*//
 
 
  return 0 ;
}
 
 
double f( const int j, const double x)
{
  double tmp ;
 
  switch(j){
  case 0 :
    tmp = 1.0 ;
    break ;
  case 1 :
    tmp = cos(x) ;
    break ;
  case 2 :
    tmp = pow( x, 3.0) ;
    break ;
  default :
    tmp = 0.0 ;
    break ;
  }
 
  return tmp ;
}

Result

polynomial_fit.png