## 第8章 非線形パラメータ推定システム

## 1. NPE の概要

NPE (Nonlinear Parameter Estimation module) は，最適化法を利用して最適化問題の解を求めるためのシステムです．NPE は，実験データからモデルパラメータを推定したり，要求された条件を満たすフィルタ，制御系，電気回路を設計するなど，モデリングにおいて大きな力を発揮します．また，NPE では，最適化法を利用して問題を解く際に必要となる付随的な作業を極力軽減し，最適化法のプログラミングから推定実行までを全て自動化しています．NPE の特徴としては，

- 統計学，数値解析の知識は特に必要なく，最適化法の特徴を知るだけで容易に利用できる．

- 非制約，制約条件付き最適化問題を扱うことができる．

- Simplex 法，BFGS 法，DFP 法，SSVN 法，共役勾配法といった代表的な最適化法を装備しており，用途に応じて選択することができる．

- プログラミング言語 (現在は，C 言語，C++ 言語をサポート) で記述されたモデル以外にも，SATELLITE に統合されている神経回路シミュレータ NCS を利用して記述されたモデルも扱うことができる．

などが挙げられます．

C 言語によるモデルを用いたパラメータ推定

NPE を利用したパラメータ推定は，(1) モデルの記述 (項2. 「モデル記述」)，(2) フィッティングに用いるデータファイルの準備 (項3. 「データファイルの準備」)，(3) 推定条件の設定 (項4. 「推定条件の設定」)，(4) 推定実行 (項5. 「NPE の実行」)，という手順で行われます．以下では，こうした手順に従ったパラメータの推定法を，

という 4 つのパラメータ ( ) をもつ 2 出力 ( ) モデルを例に説明していきます．

## 2. モデル記述

NPE を用いてパラメータ推定を行うためには，対象となるモデルを NPE の仕様に従って記述しておく必要があります．また，制約条件付きの推定を行う場合には，制約条件を設定するプログラムも決められた入出力関係を満たすように記述しなければなりません．本節では，NPE で扱うことのできるモデルと制約条件の記述形式について説明します．

### 2.1. C 言語によるモデル記述

NPE では，C 言語もしくは C++ 言語を用いてモデルを記述することができます．モデル関数では，時刻に対するデータインデックスおよびパラメータを受け取り，これらに基づいてモデル出力を計算するように記述します．このとき，モデル関数の引数は以下のように対応づけられます．

第 1 引数 : 時刻に対応するデータインデックス (int 型の変数)
第 2 引数 : 推定したいパラメータ (double 型の配列変数)
第 3 引数 : モデル出力 (double 型の配列変数)

第 2 引数のパラメータは，項4. 「推定条件の設定」 で説明する初期値設定関数 cinit() で宣言した順番で引き渡されますので，モデル記述の際はその対応に注意してください．また，第 3 引数のモデル出力もフィッティングに用いるデータファイルに格納されている順番に対応しているため，複数の出力を持つモデルを扱う場合には，配列変数に代入するモデル出力の順番に注意が必要です．

この形式に従って式 (8.1)，(8.2) の例を記述したモデルファイルをリスト 8 に示します．例のように，モデル関数内で時刻が必要な場合は，リスト 8 の 14 行目にあるように get_sampling() 関数で SATELLITE で設定されているサンプリング周波数を利用できますので，point/samp を計算することによりデータインデックス point に対応する時刻をモデル関数内で利用することができます．

ここでは例として単純な関数を取り上げましたが，記憶・内部状態を持つ連続系モデルを扱いたい場合は，モデル関数内部で static 変数などを利用するか，次節で説明する NCS 言語を用いて記述してください．

>     1   /*  NPE USRモデルファイル (sin2.cpp)  */
>     2   /*  SINの振幅, 位相, バイアス を求める  */
>     3   /*  COSの振幅, 位相, バイアス を求める  */
>     4   /*  2-出力の場合  */

>     5   #include <stdio.h>
>     6   #include <math.h>
>     7   #include <libsatellite.h>
>     8   #if !defined(PI)
>     9   #define PI 3.14159
>     10  #endif

>     11  int model(int point, double *param, double *output)
>     12  {
>     13      double y, theta, samp;
>     14      samp = get_sampling();
>     15      theta = 2.0 * PI * ((double)point / samp) + param[2];
>     16      y = param[0] * sin(theta) + param[3];
>     17      output[0] = y;
>     18      y = param[1] * cos(theta) + param[3];
>     19      output[1] = y;
>     20      return(0);
>     21  }

リスト 8 : C 言語によるモデルの記述例 (sin2.cpp)

### 2.2. NCS 言語によるモデル記述

NPE では，NCS 言語で記述されたモデルに対してもパラメータを推定できます．但し，この場合も C 言語でモデル記述した場合と同様に，モデルファイルをコンパイルしてオブジェクトファイルにしておく必要があります．NCS モデル (例 : sample.mdl) は，以下のコマンドを実行することによりコンパイルできます．

>     []SATELLITE[]~/home/demo:[1]% npp ("sample")
>     []SATELLITE[]~/home/demo:[2]% nlink ("-O2","O")

### 2.3. 制約条件の記述

NPE では，推定したいパラメータの取り得る範囲があらかじめわかっている場合，そうした情報をペナルティ関数として C 言語で記述しておくことで，制約条件付きの推定ができます．式 (8.1)，(8.2) の例において，各パラメータの範囲を

と定めた場合，制約条件の記述はリスト 9 のようになります．ペナルティ関数では，C 言語で記述したモデル関数の第 2 引数と同様に，初期値設定関数 cinit() で宣言した順番でパラメータ値が格納された配列変数を引数として受け取ります (リスト 9 の 6 行目の param)．また，制約条件に基づいて計算されたペナルティ (double 型) が関数の戻り値となります．この戻り値は，評価関数に加算されるため，結果的にパラメータ値が制約条件に反しないように推定が進行します．作成された制約条件プログラムは，モデルと同様，コンパイルしてオブジェクトファイルにしておく必要があります．

>     1   /*  2-出力モデル  制約条件プログラム (sin2_pena.cpp)  */

>     2   #include <stdio.h>
>     3   #include <math.h>
>     4   #ifdef __cplusplus

>     5   #define  Pena_Val   1000.0

>     6   double penalty( double param )
>     7   {
>     8       double pena = 0.0;
>     9       if( param[0]<=5. )        pena += Pena_Val;
>     10      if( param[1]<=5. )        pena += Pena_Val;
>     11      if( param[2]=<M_PI/12. )  pena += Pena_Val;
>     12      if( param[2]<=M_PI/2. )   pena += Pena_Val;
>     13      if( param[3]<=0 )         pena += Pena_Val;
>     14      return(pena);
>     15  }

リスト 9 : 制約条件の記述例 (sin2_pena.cpp)

## 3. データファイルの準備

NPE では，実験などで得られたデータとモデル出力との誤差で表される評価関数が最小になるようにパラメータを推定します．そのためには，推定を実行する前に，フィッティングに用いるデータ (以下，推定用データと呼ぶ) を格納したファイルを用意しておく必要があります．また，評価関数に任意の重み付けをする場合は，重みデータを格納したファイルも用意する必要があります．両者とも SATELLITE のファイル形式で，モデル出力と同じレコード数のデータとして用意します．ファイルには，レコード 0 からデータを格納します．例えば，モデルの出力が 2 の場合には，

>     []SATELLITE[]~/home/demo:[3]% data1 = $"data.dat":[0]
>     []SATELLITE[]~/home/demo:[4]% data2 = $"data.dat":[1]

という操作により，出力 1，2 のデータが，それぞれ，data1，data2 に取り出せるような形式のデータファイルとして用意しておく必要があります．

ここでは，リスト 8 に示した 2 出力のモデルにおいて，各パラメータの真値を と定め，適当な雑音を重畳させたものを実験データの代わりとして SATELLITE 上で作成し，ファイルとして保存する手続きを説明します．データポイントは 1000 点，サンプリング周波数は 1kHz とします．すなわち，1 秒間のデータファイルが作成されることになります．重みデータは，簡単なものにするため，単純に全て 1 としておきます．リスト 10 に，これらのデータファイルを作成するためのバッチファイルを示します．このバッチにより作成されたデータは，図 8.1. 「推定用データ例」 のようになります．黒色が式 (8.1) ならびにリスト 8 の 16 行目に記述したサイン関数，灰色が式 (8.2) ならびにリスト 8 の 18 行目に記述したコサイン関数につての推定用データです．バッチファイル内における各 SATELLITE コマンドの使用方法に関しては，リファレンスマニュアルならびにユーザーズマニュアルの該当項目を参照してください．

>     1   const dpt = 1000;
>     2   const A1  = 2;
>     3   const A2  = 3;
>     4   const B_  = PI/8;
>     5   const C_  = 2;

>     6   series y[dpt],wgt[dpt];

>     7   # sin,cos
>     8   t = (0~(dpt - 1)) / dpt;
>     9   y:[0] = A1 * sin(2 * PI * t + B_) + C_;
>     10  y:[1] = A2 * cos(2 * PI * t + B_) + C_;

>     11  # d
>     12  ran0 = nrand(dpt,1,0,0.01);
>     13  ran1 = nrand(dpt,3,0,0.01);

>     14  # f
>     15  y:[0] = y:[0] + ran0;
>     16  y:[1] = y:[1] + ran1;

>     17  #
>     18  wgt:[0] = (0~(dpt - 1)) * 0 + 1;
>     19  wgt:[1] = wgt:[0];

>     20  #
>     21  $"sin2.dat" = y;
>     22  $"sin2.wgt" = wgt;

リスト10 : データファイルの作成用バッチファイル (sin2_mkdata.sl)

図 8.1. 推定用データ例

以上の例では，推定用データを SATELLITE で作成しましたが，実際には，各々の実験環境において測定されたデータのフォーマットを，SATELLITE データに変換する処理が必要になります．各データ点が，スペース，あるいは，タブ，カンマ，改行で区切られ，アスキー形式で保存されているデータは，SATELLITE のユーティリティモジュールに登録されている text2buffer() を用いることにより，簡単に SATELLITE データに変換できます．例えば，リスト 8 に示したモデルのパラメータ推定に利用したい実験データが，

>     2.845, 4.746, 2.875, 4.845, 2.751, 4.808, 2.747, 4.779 ...

のように，ファイル "data.txt" に合計 2000 点格納されており，奇数番目のデータをモデル出力 1，偶数番目のデータをモデル出力 2 の推定用データとして利用したい場合には，

>     []SATELLITE[]~/home/demo:[5]% data = text2buffer("data.txt")
>     []SATELLITE[]~/home/demo:[6]% data
>     [0]:%   2.845   4.746   2.875   4.845   2.751
>     [5]:%   4.808   2.747   4.779   ...
>     []SATELLITE[]~/home/demo:[7]% data = trans(reform(data,(1000,2)))
>     []SATELLITE[]~/home/demo:[8]% data:[0]
>     [0]:%   2.845   2.875   2.751   2.747   ...
>     []SATELLITE[]~/home/demo:[9]% data:[1]
>     [0]:%   4.746   4.845   4.808   4.779   ...
>     []SATELLITE[]~/home/demo:[10]% $"data.dat" = data

とすれば，推定に必要なデータファイルを作成できます (章 9. データ変換システム)．ただし，実験データがバイナリ形式でファイルに保存されている場合には，そのフォーマットを解読した上で SATELLITE データに変換するプログラムを独自に作成しなければなりません．

## 4. 推定条件の設定

NPE では，モデルのパラメータ推定を行う前に，表 8.1. 「推定条件の設定コマンド」 に示すような各種の推定条件を設定する必要があります．このうち，(1) ～ (9) は必ず設定しなければならないものです．

表 8.1. 推定条件の設定コマンド

推定条件 | 設定コマンド
--- | ---
(1) 最適化アルゴリズム | cmethod()，clsearch()
(2) 推定モデル | cmodel()
(3) 推定の制約条件 | cpenalty()
(4) 推定パラメータ | cinit()
(5) スケーリング法 | cscale()
(6) 推定の停止条件 | clogic()，cterm()
(7) 推定モデルの出力数 | cnumber()
(8) 推定用データ | cpoint()，cdata()，cweight()
(9) 推定結果の格納 | creault()，chistory()
(10) 表示する評価関数値の計算方法 | cdisp()
(11) 評価関数の計算方法 | cnorm()

表 8.1. 「推定条件の設定コマンド」 の右側は，NPE で各推定条件を設定するための関数一覧です．各関数の引数の説明は，リファレンスマニュアルを参照してください．以下では，最低限設定しなければならない推定条件について解説します．

1. 最適化アルゴリズム…cmethod()，clsearch()

   非線形最適化法に用いる最適化アルゴリズムを指定します．NPE には，最適化アルゴリズムとして，Simplex 法，BFGS (Broyden-Fletcher-Goldfarb-Shnno) 法，DFP (Davideon-Fletcher-Powell) 法，SSVM (Self-Scaling Variables Metric) 法，共役勾配法 (Fletcher-Reeves の更新式) ，共役勾配法 (Polak-Ribiere-Polyak の更新式) の 6 種類が用意されています．Simplex 法以外の最適化アルゴリズムでは，直線上で評価関数値を最小にする一次元探索法を指定する必要があります．一次元探索法には，黄金分割法と三次補間法が用意されています．

2. 推定モデル…cmodel()

   推定モデルを指定します．モデルが C 言語で記述されている場合には，第 2 引数にモデルのファイル名を指定します．NCS 言語で記述されたモデルの場合には，オブジェクトファイル名を指定します．

3. 推定の制約条件…cpenalty()

   推定の制約条件を記述したペナルティ関数のオブジェクトファイルを指定します．

4. 推定パラメータ…cinit()

   推定するパラメータの初期値などを指定します．

5. スケーリング法…cscale()

   スケーリング法を指定します．推定すべきパラメータのオーダーが不揃いの場合，収束するまでのアルゴリズムや一次元探索の反復回数が増加し，推定時間が増大していまいます．こうした問題に対して，スケーリングでは，推定前に書くパラメータに適当な係数を乗じてパラメータのオーダーを揃える操作を行い，推定速度の向上を図ります．推定終了後には，逆の操作によってパラメータを元のオーダーに戻します．

6. 推定の停止条件…clogic()，cterm()

   推定の停止条件を指定します．停止条件には，最大反復回数，評価関数値，Simplex 法の評価関数値の標準偏差が設定できます．また，clogic() は，必ず cterm() の前に実行しておく必要があります．

7. 推定モデルの出力数…cnumber()

   推定モデルの出力数を指定します．この関数で指定した数は，推定用データや重みデータを格納しているファイルのレコード数と一致している必要があります．

8. 推定用データ…cdata()，cweight()，cpoint()

   推定に用いるデータファイルを指定します．cdata()，cweight() では，それぞれ推定用データ，重みデータを格納したファイル名を指定します．また，cpoint() ではデータファイル中のデータ点数を指定します．

9. 推定結果の格納…cresult()，chistory()

   推定結果を格納するファイルを指定します．cresult() ではモデル出力の履歴を，chistory() でパラメータ・評価関数値の履歴を格納するファイルをそれぞれ指定します．

このような推定条件をリスト 8 のモデルに対して設定する SATELLITE のバッチファイルの例 (sin2_com.sl) をリスト 11 に示します．ここでは，最適化法として経験的に優れているといわれている BFGS 法を指定し，一次元探索法には黄金分割法を選択しています．また，スケーリングは無しに設定しています．

こうしたバッチファイルは SATELLITE のシェルコマンドに登録されている inline() により，

>     []SATELLITE[]tom//home/demo:[11]% inline("sin2_com.sl")

とすることで実行できます．また，一度設定した推定条件は，関数 cstore() を用いることで

>     []SATELLITE[]tom//home/demo:[12]% cstore("sin2.prm")

のようにファイル (sin2.prm) に格納することもできます．このパラメータファイルを関数 cload() により読み込んで推定条件を同様に設定することができます．

>     1   ## NPE 推定条件スクリプト(sin2_com.sl) (USRモデル)

>     2   # 最適化アルゴリズム
>     3   cmethod("bfgs");               # "BFGS法"を設定
>     4   clsearch("golden",0.00001);    # 直線探索法として"黄金分割法"を設定

>     5   # 推定モデル
>     6   cmodel("USR","sin2.cpp");      # C言語のモデル"sin2.cpp"を設定

>     7   # 推定の制約条件
>     8   cpenalty("sin2_pena.cpp");     # ペナルティ関数"sin2_pena.cpp"を設定

>     9   # 推定パラメータ
>     10  cinit(4,1,1.8,"VAR","A1",0.000001); # "A1"の初期値を"1.8"に設定
>     11  cinit(4,2,2.8,"VAR","A2",0.000001); # "A2"の初期値を"2.8"に設定
>     12  cinit(4,3,-0.1,"VAR","B",0.000001); # "B"の初期値を"-0.1"に設定
>     13  cinit(4,4,3.0,"VAR","C",0.000001);  # "C"の初期値を"3.0"に設定

>     14  # スケーリング法
>     15  cscale(0);                 # "スケーリングなし"に設定

>     16  # 推定の停止条件
>     17  clogic("0|1");             # 条件に"最大反復回数または評価関数値"を設定
>     18  cterm(0,100);              # 最大反復回数 ≧ 100
>     19  cterm(1,20);               # 評価関数値 ≦ 20

>     20  # 推定モデルの出力数
>     21  cnumber(2);                # 出力数"2"に設定

>     22  # 推定データ
>     23  cpoint(1000);              # データ点数を"1000"に設定
>     24  cdata("sin2.dat");         # 実験データ"sin2.dat"を設定
>     25  cweight("sin2.wgt");       # 評価関数の重みデータ"sin2.wgt"を設定

>     26  # 推定結果の格納
>     27  cresult("sin2_res.dat");       # モデル出力の履歴を格納するファイルを

>     28  # "sin2_res.dat"に設定
>     29  chistory("sin2_hist.dat",10);  # パラメータ・評価関数値の履歴を格納する
>     30                                 # ファイルを"sin2_hist.dat"に設定

>     31  ## ※以下のコマンドは必要に応じて設定する

>     32  # 表示する評価関数値の計算方法
>     33  cdisp(0);                  # 表示する評価関数値の計算を"重みなし"に設定

>     34  # 評価関数の計算方法
>     35  cnorm(2);                  # 評価関数の計算方法を"2-ノルム"に設定

リスト 11 : 推定条件の設定バッチファイル例 (sin2_com.sl)

## 5. NPE の実行

ここまでの説明で，パラメータ推定を実行するために必要な，モデル (sin2.cpp)，制約条件 (sin2_pena.cpp)，推定用 (sin2.dat) ならびに重みデータ (sin2.wgt) を格納したファイルの作成 (sin2_mkdata.sl) や，推定条件の設定 (sin2_com.sl) が完了しました．こうした準備が間違いなくできていれば，以下のように推定実行関数 npe() を実行するだけで，自動的にパラメータ推定が行われます (sin2_est.sl) ．

>     []SATELLITE[]~/home/demo:[13]% series result,hist
>     []SATELLITE[]~/home/demo:[14]% npe(result,hist)
>     Making Executable file is done.

>     [   NPE]   BFGS

>           <<< 0   Step:   Error = 3495.14  Penalty = 1000  >>>
>               A1          A2           B             C
>               1.8         2.8          -0.1          3

>           ....................................................

>           <<< 1   Step:   Error = 799.14   Penalty = 0  >>>
>               A1          A2           B             C
>           1.79589     2.78321      0.52453      2.55808

関数 npe() の第 1 引数 (result) には，最終的に推定されたパラメータでのモデル出力，第 2 引数 (hist) には，推定されたパラメータと誤差の履歴が格納されます．これらの内容は，推定条件設定関数 cresult()，chistory() で指定したファイルに格納されているデータと同じです．パラメータの推定状況をみると， には制約条件を超える処理値 (-0.1) を設定したので，初期状態で "Penalty=1000" とペナルティが課せられていることがわかります．しかし，1 回目の推定後には，制約条件内に全てのパラメータが更新され，ペナルティが取り消されています．

このパラメータ推定は，

>           <<< 11  Step:   Error = 799.14   Penalty = 0  >>>
>               A1          A2           B             C
>          1.97965     2.97731      0.39562       1.9979

>       Error < 20
>           3 points are stored in "sin2_hist.dat"
>       Estimation is done.

と，推定条件で設定したように，誤差が 20 以下になると 11 回の更新で終了します．パラメータの推定後に，

>     []SATELLITE[]~/home/demo:[15]% hist
>     [0]:[0]%   1.8     2.8     -0.1     3       3495
>     [1]:[0]%   1.969   2.967   0.3864   1.997   20.44

と入力すると，パラメータの履歴が確認できます．最後の行に表示されているのが最終的に推定されたパラメータです．左から， ，誤差の順に表示されています．推定用データを生成するときに用いた値は，それぞれ， でしたので，真の値に近いパラメータが推定されていることがわかります．図 8.2. 「推定結果例」 に，推定されたパラメータを用いたモデル出力 (白色) と，推定用データ (黒色) の重ね書きを示します．パラメータ推定後のモデル出力が，推定用データとほとんど重なっており，パラメータが良好に推定されていることがわかります．

図 8.2. 推定結果例

## 6. NCS 言語を用いたパラメータ推定

NPE では，NCS 言語で記述されてモデルに対してもパラメータを推定できます．以下，Hodgkin-Huxley 型モデルによるスパイク応答波形から，モデルのナトリウム電流コンダクタンス ( ) と膜容量 ( ) を推定するために NPE を適用した例をしめします．なお，NCS 言語に関しては 章 7. 神経回路シミュレータ を参考にしてください．

- モデル記述

  NCS モデルには 7 章のリスト 3 : NCS 言語による Hodgkin-Huxley モデルを用います．このモデルを保存したファイルを hh.mdl とします．

- データ作成

  実践的には実験データを NPE で扱えるデータフォーマットに変換して使用しますが，ここでは，NCS で計算したシミュレーション結果を NPE の教師データとして利用することにします (リスト) ．NCS で HH モデルをシミュレーションし，結果をファイルオブジェクトに保存します．重みファイルも作成します．また，2-3 行目で NCS モデルをコンパイルしていますが，3 行目の第 2 引数に注意してください．NCS モデルを NPE で使用する場合には，C 言語モデルとは異なり，モデルのオブジェクトファイルを用意する必要があります．nlink 関数の第 2 引数に "O" オプションを渡すことにより，オブジェクトファイルの生成を行います．

  >     1   nassign("hh");
  >     2   npp();
  >     3   nlink("-O2","O");
  >     4   ntime(10,0.01,0.01,1);
  >     5   nstim("HH",0,"P",1,0,100,3,999);
  >     6   series V;
  >     7   nout(V,"HH",0,1);
  >     8   ninteg("F");
  >     9   ncal();
  >     10  datapt = length(V);
  >     11  series vd[datapt];
  >     12  vd:[0] = V;
  >     13  wgt:[0] = (1~(length(V))) * 0 + 1;
  >     14  $"hh.dat" = vd;
  >     15  $"hh.wgt" = wgt;

  リスト 12 : 教師データ作成スクリプト (NCS によるシミュレーション : hh_mkdata.sl)

- ペナルティ関数

  C 言語のモデルのときと同じように，推定するパラメータに対する制約条件を設定することができます．ここでは C++ 言語で制約条件を設定することにします．推定するパラメータである GNa は 0 から 1000 の間，Cm は 0 から 10 の間であるという制約条件の場合には次のリスト 13 のようなファイルを作成します．

  >     1   #include <stdio.h>
  >     2   #include <stdlib.h>
  >     3   #include <math.h>
  >     4   #define lambda  0.1

  >     5   #ifdef __cplusplus
  >     6   extern "C" {
  >     7   #endif

  >     8   double penalty(double *param)
  >     9   {
  >     10     double ret=0.0;
  >     11     if(param[0] < 0. ) ret += param[0]*param[0];
  >     12     if(param[0] < 1000.) ret += (param[0]-1000.)*(param[0]-1000.);
  >     13     if(param[1] < 0. ) ret += param[1]*param[1];
  >     14     if(param[1] < 10.) ret += (param[0]-10.)*(param[0]-10.);

  >     15     return(ret * lambda);
  >     16  }
  >     17  #ifdef __cplusplus
  >     18  }
  >     19  #endif

- NPE 条件設定

  C 言語のモデルのときと同じように，パラメータを推定するときの NPE の設定条件をリスト 14 に示します．

  >     1   cmethod("bfgs");
  >     2   clsearch("golden",0.00001);
  >     3   cmodel("NCS","hh.o");
  >     4   cpenalty("hh_pena.cpp");
  >     5   cinit(2,1,140.0,"VAR","GNa@HH",0.000001);
  >     6   cinit(2,2,0.1,"VAR","Cm@HH",0.000001);
  >     7   cscale(3);
  >     8   clogic("0|1");
  >     9   cterm(0,20);
  >     10  cterm(1,0.0001);
  >     11  cnumber(1);
  >     12  cdata("hh.dat");
  >     13  cweight("hh.wgt");
  >     14  cresult("hh_res.dat");
  >     15  chistory("hh_hist.dat",10);
  >     16  cdisp(0);
  >     17  cnorm(2);

  リスト : 14 Hodgkin-Huxley モデルパラメータ推定のための設定条件 (hh_com.sl)

  C 言語を用いた場合との相違点をまとめると以下のようになります．

  - モデル指定

    3 行目
    cmodel ("NCS","hh.o")
    第 1 引数は，C 言語のモデルの場合には "USER" でしたが，NCS 言語の場合には "NCS" になります．また，第 2 引数には nlink 関数で作成したオブジェクトファイルを指定します．

  - パラメータの初期値設定

    5，6 行目
    cinit(2, 1, 140.0, "VAR", "GNa@HH", 0.000001);
    cinit(2, 2, 0.1, "VAR", "Cm@HH", 0.000001);
    推定するパラメータを指定する cinit の 5 番目の引数の書式は NCS のパラメータ名 @NCS のモジュール名となります．

- 推定の実行

  ここまで推定に必要なモデル，データ，NPE 条件ファイルがそろいましたので，あとは NPE を実行するだけです．推定のためのスクリプトをリスト 15 に示します．NCS モデルを使ったパラメータ推定の場合，NPE が連続系のシミュレーションを NCS と協調して行うわけですが，そのため，npe() 関数を実行する前に，シミュレーション条件を推定したいデータを取得した条件に合わせて設定する必要があります (リスト 15 : 1 から 8 行目)．

  >     1   nassign("hh");
  >     2   npp();
  >     3   nlink("-O2");
  >     4   ntime(10,0.01,0.01,1);
  >     5   nstim("HH",0,"P",1,0,100,3,999);
  >     6   series V;
  >     7   nout(V,"HH",0,1);
  >     8   ninteg("F");
  >     9   inline("hh_com.sl");
  >     10  cpoint(length($"hh.dat":[0]);
  >     11  series res,hi;
  >     12  npe(res,hi);

  リスト 15 : 推定実行スクリプト (hh_est.sl)

Last updated: 2005/11/12