GSLを使ってフィッティングしてみた

http://gonzaburou.cocolog-nifty.com/blog/2006/11/gslvisual_cc_8716.html
前述のとおりGSLをwindowsで使えるようにした。
あとは実際にどんなソースをかけばいいのかであるが、これは公式リファレンスの日本語訳
http://www.cbrc.jp/~tominaga/translations/#gsl
のサンプルプログラムを参考にするのがよいであろう。
私の場合は非線形二次近似をしたかったのでその項を参照。

今回試しにAl（励起状態）->Mgの半減期測定したデータの解析をしてみた。
MCS001.csv
に実験データが収められている。

/**
 * @file GSLで非線形二次近似のサンプル
 *
 * @author Hoshimi's Works
 * @date 2008/09/21
 * @note ご自由にお使いください。ただし責任は一切負いません。負えません。
 */
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <vector>
//GSL
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
//GSL
#pragma comment(lib,"../lib/gsl/gsl.lib")
#pragma comment(lib,"../lib/gsl/gslML.lib")

#define FILEPATH "./MCS001.csv"	//読み込むファイルのパス

//半減期のFitting
//求めるパラメータは
//f(x) = Aexp(-lambda t)

//この形式は固定
struct data{
	size_t n;
	double * x;
	double * y;
	double sigma;		//標準偏差
};

/**
 * @brief 更新したパラメータに対する理論値と実験値のズレの計算
 */
int exp_f(const gsl_vector * vec,void * params,gsl_vector* f){
	data* pData = (data*)params;
	size_t n = pData->n;
	double* x = pData->x;
	double* y = pData->y;
	double sigma = pData->sigma;

	double A =gsl_vector_get(vec,0);
	double lambda = gsl_vector_get(vec,1);

	for(int i = 0;i < n;i++){
		double t = x[i];							//横軸
		double Yi=A * exp(-lambda * t);				//縦軸（理論値）
		gsl_vector_set(f,i,(Yi - y[i])/sigma);
	}

	return GSL_SUCCESS;
}
/**
 * @brief 更新したパラメータに対するヤコビアンの計算
 */
int exp_df(const gsl_vector * vec,void * params,gsl_matrix* J){
	//パラメータをキャスト
	data* pData = (data*)params;

	size_t n = pData->n;
	double* x = pData->x;
	//double* y = pData->y;
	double sigma = pData->sigma;

	double A =gsl_vector_get(vec,0);
	double lambda = gsl_vector_get(vec,1);

	for(int i = 0;i < n;i++){
		double t = x[i];					//横軸
		double s = sigma;					//標準偏差
		double e = exp(-lambda * t);				//exp

		double dfdA = e/s;					//\frac{dA}{dy}/sigma
		double dfdLambda = -A*t*e/s;				//\frac{df}{d\lambda}/sigma
		//Jacobianの設定
		gsl_matrix_set(J,i,0,dfdA);
		gsl_matrix_set(J,i,1,dfdLambda);
	}
	return GSL_SUCCESS;
}

int exp_fdf(const gsl_vector * x,void * params,gsl_vector* f,gsl_matrix* J){
	return exp_f(x,params,f) | exp_df(x,params,J);
}

int main(int argc, char** argv)
{
	//外部ファイルからデータの読み込み
	FILE* fp = fopen(FILEPATH,"r");

	char line[ 1024 ];
	//char cx[8],cy[8];
	int ix=0,iy=0,nSum=0;
	using namespace std;
	vector<vector<double>> vValue;
	while( fgets( line, 1024, fp ) != NULL ){
		sscanf(line,"%d,%d",&ix,&iy);
		vector<double> vec;
		vec.push_back(ix);
		vec.push_back(iy);
		nSum += ix + iy;
		vValue.push_back(vec);
	}

	//fittingの準備
	const gsl_multifit_fdfsolver_type *T;
	gsl_multifit_fdfsolver *s;

	const size_t n = vValue.size();	//サンプルの個数
	const size_t p = 2;				//パラメータの個数

	//標準偏差の計算
	double nAve = nSum / n;
	long double nDeviation = 0;			//標準偏差
	for( int i = 0 ; i < n ;i++){
		nDeviation += pow(vValue[i][1] - nAve,2);
	}
	nDeviation /= n;

	double* pX = new double[n];
	double* pY = new double[n];			//実験値
	double sigma = sqrt(nDeviation);	//エラーバーの幅

	for( int i = 0 ; i < n ;i++){
		pX[i] = vValue[i][0];
		pY[i] = vValue[i][1];
	}

	data d = {n,pX,pY,sigma};

	gsl_multifit_function_fdf f;

	double x_init[2] = {80000.0,10.0};	//パラメータの初期化

	gsl_vector_view x = gsl_vector_view_array(x_init,p);

	const gsl_rng_type * type;
	gsl_rng * r;

	gsl_rng_env_setup();
	type = gsl_rng_default;
	r = gsl_rng_alloc(type);

	f.f = &exp_f;
	f.df = &exp_df;
	f.fdf = &exp_fdf;
	f.n = n;
	f.p = p;
	f.params = &d;

	T = gsl_multifit_fdfsolver_lmsder;
	s = gsl_multifit_fdfsolver_alloc(T,n,p);
	gsl_multifit_fdfsolver_set(s,&f,&x.vector);

#define GETV(i) gsl_vector_get(s->x,i)

	int status = GSL_CONTINUE;
	for(int counter = 0;counter < 500 && status == GSL_CONTINUE;counter++){
		status = gsl_multifit_fdfsolver_iterate(s);

		printf("%f %f\t%f\n",GETV(0),GETV(1),gsl_blas_dnrm2(s->f));

		if(status){break;}

		status = gsl_multifit_test_delta(s->dx,s->x,1e-4,1e-4);
	}

	gsl_matrix * covar = gsl_matrix_alloc(p,p);	//行列の確保
	gsl_multifit_covar(s->J,0.0,covar);

#define FIT(i) gsl_vector_get(s->x,i)
#define ERROR(i) sqrt(gsl_matrix_get(covar,i,i))


	printf("A\t=%.5f +/- %.5f\n",FIT(0),ERROR(0));
	printf("lambda\t=%.5f +/- %.5f\n",FIT(1),ERROR(1));

	double chi = gsl_blas_dnrm2(s->f);
	printf("chisq/dof = %d",pow(chi,2.0)/(n-p));

	delete[] pX;
	delete[] pY;

	int end;
	std::cin>>end;

	return 0;
}

出力結果は以下のとおり

80004.462769 6.916877	54.158680
80013.396685 6.634422	54.158598
...(略)...
115798.297218 0.002936	5.429567
115796.485374 0.002936	5.429567
A	=115796.48537 +/- 3317.91143
lambda	=0.00294 +/- 0.00015
chisq/dof = 321559052

グラフに出力してみた

…これはフィッティングできてるのかな？かな？

に対する近似にしときゃよかったと後で後悔。いや、これからやりなおしますけどね。
…できていることにしよう。最初のほうはなんかゴミはいってるっぽいし。