GSL関数をフィット関数として使用する方法を説明したものです。
必要なOriginのバージョン:8.0 SR6
1.次の関数で、下にあるサンプルデータをフィットします。
0.1 0.10517 0.2 0.2214 0.3 0.34986 0.4 0.49182 0.5 0.64872 0.6 0.82212 0.7 1.01375 0.8 1.22554 0.9 1.4596 1 1.71828 1.1 2.00417 1.2 2.32012 1.3 2.6693 1.4 3.0552 1.5 3.48169 1.6 3.95303 1.7 4.47395 1.8 5.04965 1.9 5.68589 2 6.38906 2.1 7.16617 2.2 8.02501 2.3 8.97418 2.4 10.02318 2.5 11.18249 2.6 12.46374 2.7 13.87973 2.8 15.44465 2.9 17.17415 3 19.08554 3.1 21.19795 3.2 23.53253
2.次のステップに進む前に、Originをインストールしたフォルダのユーザファイルフォルダにocgsl.h ファイルを追加します。同じ場所に、Calling GNU Scientific Libraryからgsl_dllファイルをコピーしてください。
ocgsl.h
#pragma dll(libgsl, header) // これはOCの特別なプラグマ で、 // ヘッダキーワードはlibgsl.dllがこのファイルと同じ場所にあることを示しています。 #define GSL_EXPORT // OCでは、これは不要ですので、削除します。 //gsl関数のプロトタイプをここで直接探すことができます。 typedef double (* FUNC)(double x, void * params); struct gsl_function_struct { FUNC function; void * params; }; typedef struct gsl_function_struct gsl_function ; typedef struct { size_t limit; size_t size; size_t nrmax; size_t i; size_t maximum_level; double *alist; double *blist; double *rlist; double *elist; size_t *order; size_t *level; } gsl_integration_workspace; GSL_EXPORT gsl_integration_workspace *gsl_integration_workspace_alloc (const size_t n); GSL_EXPORT void gsl_integration_workspace_free (gsl_integration_workspace * w); GSL_EXPORT int gsl_integration_qag (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, int key, gsl_integration_workspace * workspace, double *result, double *abserr);
3.F9を押し、フィット関数オーガナイザを開き、以下のように新しい関数を定義します。
4.関数フィールドの右側にあるボタンをクリックし、コードビルダを開き、以下のコードを追加して、 _nlfgsl_integration_qag.fitをコンパイルします。
#include "..\ocgsl.h" static double f_callback(double x, void * params) { double alpha = *(double *)params; return exp(alpha*x); } void _nlsfgsl_integration_qag( // フィットパラメータ double y0, double a, double beta, // 独立変数: double x, // 従属変数: double& y) { // 編集可能部分開始 double result, err, expected = -4.0; // 1000個の倍精度の間隔を持つことができるワークスペースを確保します。 // それらの積分結果とエラーを推定します。 gsl_integration_workspace *ww = gsl_integration_workspace_alloc(1000); gsl_function F; F.function = f_callback; F.params = &beta ; // 積分範囲 (0, x), は求められている絶対エラー０ //から、相対エラー1e-7間にあります gsl_integration_qag(&F, 0, x, 0, 1e-7, 1000, 0, ww, &result, &err); // ワークスペースwに関係したメモリが解放されます gsl_integration_workspace_free (ww); y = y0 + a*result; // 編集可能部分終了 }
さらに、以下のコードを追加すればフィット関数は完璧です。
//---------------------------------------------------------- // #include <ONLSF.h> #include "..\ocgsl.h" static double f_callback(double x, void * params) { double alpha = *(double *)params; return exp(alpha*x); } void _nlsfgsl_integration_qag( // フィットパラメータ double y0, double a, double beta, // 独立変数: double x, // 従属変数: double& y) { // 編集可能部分開始 NLFitContext *pCtxt = Project.GetNLFitContext(); if ( pCtxt ) { static vector vInteg; NLSFCURRINFO stCurrInfo; pCtxt->GetFitCurrInfo(&stCurrInfo); int nCurrentIndex = stCurrInfo.nCurrDataIndex; BOOL bIsNewParamValues = pCtxt->IsNewParamValues(); if ( bIsNewParamValues ) { vector vx; pCtxt->GetIndepData(&vx); int nSize = vx.GetSize(); vInteg.SetSize(nSize); // 1000個の倍精度の間隔を持つことができるワークスペースを確保します。 // それらの積分結果とエラーを推定します。 gsl_integration_workspace *ww = gsl_integration_workspace_alloc(1000); gsl_function F; F.function = f_callback; F.params = &beta ; double result, err, expected = -4.0; for(int ii=0; ii<nSize; ++ii) { // 積分範囲(0, vx[ii]), は求められている絶対エラー０ // から、相対エラー1e-7間にあります gsl_integration_qag(&F, 0, vx[ii], 0, 1e-7, 1000, 0, ww, &result, &err); vInteg[ii] = result; } // ワークスペースwに関係したメモリが解放されます gsl_integration_workspace_free (ww); } y = y0 + a*vInteg[nCurrentIndex]; x; } // 編集可能部分終了 }
5.次の初期化コードを追加します。
パラメータ初期化
//パラメータを初期化するコード sort( x_y_curve ); double coeff[2]; fitpoly( x_y_curve, 1, coeff); a = coeff[0]; y0 = coeff[1]; beta=1.0
6.ユーザ定義関数 gsl_integration_qagを使ってフィットすると、以下の結果が得られます。
y0 = -1.06363E-6
a = 1
beta =1