GNU Scientific Libraryを使ったユーザ定義フィット関数
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
|