|
Copyright (c) 2010 ykuroda All right reserved. |
|
Copyright (c) 2010 ykuroda All right reserved. |
１）非線形最小二乗法
SPANAで用いられている非線形最小二乗演算アルゴリズムはDamping Gauss-Newton法で、最適化するべきパラメーター・ベクトル
P の修正ベクトル ΔP に縮小因子(Damping Factor) α ( 0 < α ≦ １ ）を導入して、
P k+1 = P k + α ΔP
を求め、測定値と理論値の二乗残差(ss)を最小化するパラメータを逐次探索する。
最適解の収束判定は標準で、*
a) 二乗残差の減少量が残差の 0.001%以下
b) 各パラメータの変動値がパラメータ値の0.01％以下
になった時としている。
＊ これらの判定基準値は最適化計算設定ダイアログボックスの"Damping Factor"表示先頭にある"*"マークをクリックして現れる画面から変更することが出来る。
２）パラメータ初期値と収束
最小二乗法は二乗残差曲面の局所的な線形性を仮定して残差極小値を探索するため、真値から遠い不適切なパラメータ初期値から出発すると真のミニマムに至る前に局所ミニマムに陥って収束と判定される場合がある(図１）。そのようなケースではパラメータ初期推定値の選択が重要になる。例えば、1：1錯体生成系の滴定で適切なパラメータ初期値は、差スペクトルのプロット図などからΔAbsの外挿極限値ΔAbsmaxを推定してサンプル濃度「A]0から
Δεestimate ～ ΔAbsmax/[A]0
ΔAbsmax/2 を示す滴定試薬 B の濃度 [B]1/2 から
K estimate ～ 1/[B]1/2
として妥当な値が得られる(図２）。
それでも妥当な収束値が得られない場合の対策としては、
a) 初期推定値を大幅に変更して計算する
b) 残渣変動に対して値が過剰に反応するパラメータ（例えば冪数）を固定して計算する
c) 極端に外れたデータを除外する
ｄ）Damping factorの初期値を変更する１）
などが有効な場合がある。
図１ 
図２
３）一般データの解析
一般にSPANAの最小二乗解析では、データは読み込んだスペクトルデータから解析に必要な処理を経て"Result"ダイアログボックスに表形式で収集した数値データ (t, Absi) のみを用いて解析が行われる。したがって、この表に適当な数値データ (t, Absi等) が設定されていればSPANA由来のデータ以外のデータに対してもSPANAの最小二乗計算を適用できる。表に数値データを設定するには、"Analysis-Least
Square Analyses-User's Data" メニューから表の行・列数を設定して必要なディメンジョンの表をResultダイアログボックス上に生成させて、表に直接キーボードから数値を入力することで行うことが出来る。あるいは、データを
"*.LSQ"タイプのTextファイルにしてSPANAから直接このファイル読み込んでも良い（この"*.LSQ"ファイルのフォーマットについてはマニュアルを参照）
解析の具体例については、ＮＭＲデータの解析、濃度データの解析、などを参照。
４）ユーザー定義関数の最小二乗解析２） (Analysis, Exampleの項も参照）
SPANAには内部ルーチンとして定義されていないモデルに対しても理論値を生成する数式モデルをユーザーが自らプログラミングし、それを用いてデータを解析する機能がある。その目的で、最小二乗演算のためのルーチンとしてSPANA用にC言語で書かれた一般最小二乗最適化計算プログラム、
spana_lsq.c と、速度解析のために特化した最小二乗最適化プログラム、spana_redap.c を提供する。解析のための実行プログラムは、最小二乗計算の対象とする理論値を算出する関数（ function ）をサブルーチンとしてユーザー自身が作成し、これら最小二乗演算ルーチンとリンク・コンパイルして生成する。
例として、NMRの化学シフト変化から錯体生成定数を求めるプログラムを圧縮ファイル名"NMRTIT"としてサンプルデータ、解析例などと共に、ExternalPrograms.zip
に含めたので参照されたい（NMRTITの機能はSPANA本体のNMRデータ解析機能と実際上同じものである）。
a) 一般非線形最小二乗法 ( spana_lsq.c )
spana_lsq.c のためのサブルーチンは、ユーザーの系について測定値を従属変数（y[i]）、実験で設定された条件を独立変数（x[1]）、とするモデル関数をつくり、これを用いて理論値の数値計算を行うプログラムとして作成される。その利用にあたって、ユーザーはfunction(int
i, double x[], double p[]) として定義されるサブルーチン内に以下の予め宣言されている４つの配列変数を用いて、サブルーチン終了時に
y[i] の理論値が確定するよう演算のアルゴリズムを定義する。
<< functionサブルーチンのために予め宣言されている配列 >>
y[j] ： 従属変数 （スペクトル滴定におけるΔAbsなど測定値）、j = 1 - 20
x[1]： 独立変数 （滴定における滴定試剤濃度など"t"変数に対応し、SPANAでは配列数は"1"に固定）
c[k] ： 条件設定定数（滴定における被滴定サンプル濃度など実験中変化しない定数）、k = 1 - 30
p[l] ： 最適化すべきパラメータ、l = 1 - 20
以下に "function"サブルーチンの一般的な形式を示す。
#include <math.h>
extern double y[],c[];
void function(int i, double x[], double p[])
{
(double ......; )
(int ......; )
/* define y[i] as functions of x[1], p[j], and c[k] here */
/* i.e., y[i] = f (x[1], p[j], c[k] ) */
return;
}
最初の３行 #include <math.h> 、extern double y[],c[]、および void function(int
i, double x[], double p[]) は、このプログラムがspana_lsq.c によって呼び出されるサブルーチンであることを示すもので変更してはならない。また整数変数
i は複数の従属変数 y[]の中から個別のy[j] を指定する変数としてfunction 関数内で宣言されている。ユーザーはこれらの変数を使ってモデル関数を
function の { } 内に定義するが、上記の他に変数を使用する場合には実数は "double" 、整数は "int"
としてプログラム内で宣言する必要がある。
i) プログラム例-1
単純な例として、独立変数 x と三つのパラメータ a, b, c による2次関数で記述されるデータ値、 y = ax2 + bx + c、 に対する最適化プログラムを以下に示す。
y を y[1]、x を x[1] 、a を p[1]、b を p[2]、c を p[3]、として
#include <math.h>
extern double y[],c[];
void function(int i, double x[], double p[])
{
y[1]=p[1]*x[1]*x[1]+p[2]*x[1]+p[3];
return;
}
というテキストを作り、これを例えば "test.c"と言った"c"を拡張子とする単純テキスト・ファイルとして保存する。
これを、spana_lsq.c とリンクしてコンパイルすれば SPANA から使用できる exe 形式の実行プログラムが得られる。コンパイラは標準的なCコンパイラでよいが、例えばネット上で入手できるフリー・コンパイラであるMicrosoft Visual Studio C++ などが簡便に使用できる。
例えば Visual Studio C++ のコマンドから
cl test.c spana_lsq.c
とすると、実行プログラムtest.exeが得られ、これがSPANAから呼び出されるプログラムとなる。3)
ii) プログラム例-2
以下、次のような逐次 １：１：１ 錯体の生成系を例にとって実際のプログラム作成を考える。
A ＋ B ↔ AB K1＝[AB] / [A][B] eq.1
AB + C ↔ ABC K2= [ABC] / [AB][C] eq.2
この系で吸収をもつ A をスペクトル上透明な B と C の混合液で滴定する実験を考える。A の初期濃度を[A]0、B, C の全濃度を[B]T、[C]Tとすると滴定中の化学種の濃度関係は
[A]0 = [A] + [AB] +[ABC] eq.3
[B]T = [B] + [AB] + 「ABC] eq.4
[C]T = [C] + [ABC] eq.5
最小二乗計算の対象となる差スペクトルΔAbs は
ΔAbs = Abs - εA[A] 0 = Δε1[AB] + Δε2[ABC] eq.6
（ Δε1= εAB - εA , Δε2= εABC - εA）
実験において指定できる変数である[A]0 、[B]T 、[C]T 、及び最適化パラメータであるK1 、K2を用いて算出される AB および ABC の濃度を６式に代入すれば測定値 ΔAbs の理論値が算出される。
一般にこれら錯体の濃度は平衡式（eq.1, eq.2）と濃度関係式（eq.3 - 5）からなる連立方程式の解として得られるが、この系ではその解は２次以上の方程式の解となり簡単な解析的式として表すことが困難である。
しかし、ここで必要な解は実際の実験条件下で化学的に意味のある数値解だけなので、コンピュータによる近似数値解法を用いることが出来る。
以下、[B]の濃度を算出する手順例を示す。
５式（および１、２式）より
[C] = [C]T / ( 1 + K1K2[A][B] ) eq.7
７式（および１式）を３式に代入して、
[A]0 = [A] + K1[A][B] + K1K2[C]T[A][B] / ( 1 + K1K2[A][B] ) eq.8
が得られる。これを整理して [A] の２次方程式
K1K2[B] ( 1 + K1[B] ) [A]2 + ( 1 + K1[B] ( 1 + K2 ([C]T -[A]0 ))) [A] - [A]0 = 0 eq.9
を得る。[B] の関数としての [A] の解 [A]B は
a = K1K2[B] ( 1 + K1[B] ) , b = ( 1 + K1[B] ( 1 + K2 ([C]T -[A]0 ))) , c = -[A]0
として、
[A]B = (-b + (b2 - 4ac）1/2)/2a eq.10
より得られる。[A]B を７式に代入して[C] を[B] の関数として表した[C]B は、
[C]B = [C]T / ( 1 + K1K2 [A]B[B] ) eq.11
[A]B、[C]B を４式に代入すると
[B]T = [B] + K1[A]B[B] + K1K2[A]B[B][C]B eq.12
が得られ、[B] は１２式の解である。
未知変数として[B] のみを含む１２式を解析的に解くのは煩雑だが、K1 , K2 , [A]0 , [B]T , [C]T が正値（化学的に有意）であれば [B] は
[B]T < [A]0 の時 0 < [B] < [B]T に、
[B]T > [A]0 の時 [B]T - [A]0 < [B] < [B]T に、
一つの解を持ち、これは１２式から得られる関数、
f ([B]) = [B] + K1[A]B[B] + K1K2[A]B[B][C]B - [B]T eq.13
に対して f ([B]) = 0 を満たす[B] として得られる。
得るべき解の範囲が与えられているので、ここでは近似解を得る方法として最も簡便で応用範囲も広い２分法 (http://en.wikipedia.org/wiki/Bisection_method) を用いる。これは[B]min= 0 or [B]T - [A]0 , [B]max= [B]T から始めて、[B]mid = ([B]max + [B]min) / 2 における f ([B]mid) を評価して、
f ([B]max) x f ([B]mid) < 0 or f ([B]min) x f ([B]mid) < 0
を満たすように[B]mid を [B]max 又は [B]min と置き換えながら f ([B]mid) < ε（計算機誤差）となるまで繰り返し計算を行うことで[B] の数値解を求める方法である。
[A] および [C] は[B] の値を使ってそれぞれ eq.10 および eq.11 から、また [AB] および [ABC] はこれらの値を用いて
eq.1 および eq.2 から算出できる。
以上の手順を組み込んだサブルーチン例として、[B]T = [C]T の滴定液で滴定し、最大５点の波長でデータを収集する実験系を想定したプログラムを以下に示す。
/* Function for ABC Complexation Analysis */
/* A + B <---> AB K1= [AB]/[A][B] */
/* AB + C <---> ABC K2= [ABC]/[AB][C] */
/* Delta Abs = e1*[AB] + e2*[ABC] */
/* [A]0 = c[1], [C]T = [B]T = x[1] */
/* Delta Abs = y[i] (i = 1 - 5) */
/* K1 = p[1], K2 = p[2], */
/* e1 = p[3] - p[11], e2 = p[4] - p[12] */
#include <math.h>
extern double y[],c[];
double mfunction2(double h,double x[],double p[]);
void function(int i,double x[],double p[])
{
double max1,min1,mid,fl,fs,fm,af,bf,cf,ab,abc,aa,bb,cc;
int j;
if (x[1]==0) {ab=0;abc=0;goto org;} /* in case of [B]T = 0 */
max1=x[1]; /* start bisect algorism */
if (x[1]<c[1]) min1=0; else min1=x[1]-c[1];
fl=mfunction2(max1,x,p);fs=-x[1]; /* f([B]max) & f([B]min) */
if (fl*fs>0) goto fend1; /* goto Error treatment */
for(j=1;j<=1000;j++) /* search [B] routine */
{
mid=(max1+min1)/2;fm=mfunction2(mid,x,p); /* f([B]mid) */
if (fm*fs>=0) min1=mid;
else
{if (fm*fl>=0) max1=mid; else goto fend1;}
if (fabs(fm)<1e-35) break; /* found [B] */
}
bf=mid; /* bf (free [B]) solution */
aa=p[1]*p[2]*bf*(1+p[1]*bf); /* coefficient aa */
bb=1+p[1]*bf*(1+p[2]*(x[1]-c[1]));cc=-c[1]; /* coefficient bb & cc */
af=(-bb+sqrt(bb*bb-4*aa*cc))/2/aa; /* af (free [A]) */
cf=x[1]/(1+p[1]*p[2]*af*bf); /* cf (free [C]) */
ab=p[1]*af*bf; /* ab ([AB]) */
abc=p[2]*ab*cf; /* abc ([ABC]) */
org:
switch (i)
{
case 1:
y[1]=p[3]*ab+p[4]*abc; /* Delta Abs = e1*[AB] + e2*[ABC] */
break;
case 2:
y[2]=p[5]*ab+p[6]*abc;
break;
case 3:
y[3]=p[7]*ab+p[8]*abc;
break;
case 4:
y[4]=p[9]*ab+p[10]*abc;
break;
case 5:
y[5]=p[11]*ab+p[12]*abc;
break;
}
return;
fend1:
printf("\nFail in finding the one root of function2");
}
double mfunction2(double h,double x[],double p[])
{
double f,af,cf,aa,bb,cc;
aa=p[1]*p[2]*h*(1+p[1]*h); /* bf=h */
bb=1+p[1]*h*(1+p[2]*(x[1]-c[1]));cc=-c[1];
af=(-bb+sqrt(bb*bb-4*aa*cc))/2/aa; /* af=[A] , eq.10 */
cf=x[1]/(1+p[1]*p[2]*af*h); /* cf=[C] */
f=h+p[1]*af*h+p[1]*p[2]*af*h*cf-x[1]; /* f([B])=[B]+[AB]+[ABC]-[B]o , eq.13 */
return (f);
}
プログラム中 /* */ で囲まれた部分はコメント行で、実行プログラムとしての意味は無い。５点の測定波長におけるΔAbs はそれぞれ y[1]
- y[5]、[B]T ( = [C]T) は x[1]、最適化すべきパラメータは、K1 ⇔ p[1], K2 ⇔ p[2],、５波長でのΔε1、Δε2は p[3] - p[12] 、[A]0 は c[1] として定義されている。
また f ([B]) の評価は新たに定義した関数 double mfunction2(double h, double x[], double p[]) で行っている。
これを、spana_lsq.c とリンクしてコンパイルすれば SPANA から使用できる exe 形式の実行プログラムが得られる。コンパイラは標準的なCコンパイラでよいが、例えばネット上で入手できるフリー・コンパイラであるMicrosoft Visual Studio C++ などが簡便に使用できる。その場合は、上記ソース・プログラムをテキストファイルとして例えば"abc_comp.c"の名で保存した後、
Visual Studio C++ のコマンドから
cl abc_comp.c spana_lsq.c
とすると、実行プログラムabc_comp.exeが得られ、これがSPANAから呼び出されるプログラムとなる。3)
ｂ） 反応速度解析用最小二乗法 （ spana_redap.c ）
化学においてスペクトルを利用する測定の一つに反応速度の追跡がある。SPANAもスペクトルの経時変化から幾つかの速度過程の解析を行う機能を持っているが、ここではSPANAのユーザー定義システムに対する最小二乗解析機能を利用したより応用範囲の広い速度解析手法について説明する。
反応速度の解析において、しばしば問題になるのは反応式が確定すれば速度式は容易に書き下すことが出来るが、実際の測定データと理論式の対応に必要な経時変化の積分形理論式を得ることが困難なことがあることである。微分型の反応速度式から、積分型の理論式を得るために実験条件に制限を設けて、反応次数を下げたり、定常状態近似を適用することによる簡略化近似がしばしば有効な手法として用いられるが、こういった条件を適用することが実験的に困難であったり、あるいは、そういった条件によって情報が失われてしまう場合などもあって、常に簡略化が有効なわけではない。
SPANAで提供する速度解析用のプログラム "spana_redap.c"は微分型速度式より数値積分でデータに対応する理論値を求め最小二乗演算を行う。従って、ユーザーは反応の積分型速度式を求めることなく速度解析を行うことが出来る。以下、spana_redap.c
にリンクするユーザーのサブルーチン作成の概略を説明する。4）
ユーザーの反応系を記述するfunction 関数で使用する変数として以下の４つが予め宣言されている。
<< function サブルーチンのために予め宣言されている配列 >>
x[i] ： 反応系を構成する物理量（化学種濃度、吸光度、など）、i = 1 - 20
dx[i] ： x[i] の時間微分、dx[i]/dt
p[j] ： 最適化すべきパラメーター、j = 1 - 20
c[k]： 測定中変化しない条件設定定数、k = 1 - 10
function 関数の一般形を以下に示す。
#include <math.h>
extern double dx[],x[],c[],p[];
void function()
{
/* define dx[i]/dt */
}
spana_redap では数値積分はRunge-Kutta-Gill 法 (http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) によってメインルーチン内で行われるので、function 関数には微分方程式のみを定義すればよい。
一例として、基質 A と前段平衡で結合する触媒 C によって加速を受ける Ａ → B の反応について考える。
k+ k2
A ＋ C ↔AC → B ＋ C eq.14
k-
このタイプの反応は酵素反応のように基質濃度が触媒濃度に対して十分大きい場合は中間体 AC に対して定常状態近似が適応できてMichaelis-Menten
型の反応となることは良く知られているが、そういった条件をとることができない場合もある。
この反応の速度関係は以下のように書くことが出来る。
d[A]/dt = -k+[A][C] + k-[AC] eq.15
d[C]/dt = -k+[A][C] + (k- + k2)[AC] eq.16
d[AC]/dt = -(k- + k2)[AC] + k+[A][C] eq.17
d[B]/dt = k2[AC] eq.18
spana_redap では、例えば A : x[1]、B : x[2]、C : x[3]、AC : x[4]、k+ : p[1]、k- : p[2] 、k2 : p[3] として function関数は次のように定義される。
/* Catalytic Reaction FOR REDAP */
/* A + C <-> AC -> B + C */
/* d[B]/dt = k[AC] */
/* p[1] = k+ , p[2] = k- , p[3] = k2 */
/* x[1] = [A] , x[2] = [B] , x[3] = [C] , x[4] = [AC] */
#include <math.h>
extern float dx[],x[],c[],p[];
void function()
{
dx[1]=-p[1]*x[1]*x[3]+p[2]*x[4];
dx[2]=p[3]*x[4];
dx[3]=-p[1]*x[1]*x[3]+(p[2]+p[3])*x[4]
dx[4]=-(p[2]+p[3])*x[4]+p[1]*x[1]*x[3];
}
void function2()
{
/* Define necessary equations to determined initial x[i] */
/* which is identified as '99999'. */
/* x[i]=g(x[j],p[i],c[i]) */
}
このプログラムでは void function2() として実質 /* */ で囲まれたコメント行だけを含む空の関数が書かれているが、これは後述する平衡反応を含む反応などで変数の初期条件を算出する必要のある場合に使用する予備関数で、ここでは使用しないのでこのまま削除せず置いておいてよい。
これを例えば "cat-react.c"というテキスト・ファイルとして保存し、Microsoft Visual Studio C++ から
cl cat-react.c spana_redap.c
としてコンパイルすれば ＳＰＡＮＡ から使用することの出来る最小二乗用プログラム cat-react.exe が作られる。
この反応で前段可逆反応が全反応中で事実上平衡反応として扱うことが出来るほど後続反応に比べて十分に早い場合、
K = [AC]/[A][C] eq.19
と平衡式を置くことが出来きて反応はMichaelis-Menten 型の反応になる。この、関係式を両辺微分して
K[A]・d[C]/dt + K[C]・d[A]/dt = d[AC]/dt eq.20
また、Ａ 及び C の設定濃度（一定）をそれぞれ[A]T、[C]T として
[A]T = [A] + [AC] + [B] eq.21
[C]T = [C] + [AC] eq.22
これを両辺微分して、
d[A]/dt + d[AC]/dt + d[B]/dt =0 eq.23
d[C]/dt + d[AC]/dt = 0 eq.24
の関係が得られる。これらの関係と、生成物の生成反応プロセスに対する d[B]/dt = k2[AC] の関係より、各化学種の速度式として
d[A]/dt = -k2[AC] (1 + K[A])/(1 + K ([A] + [C])) eq.25
d[B]/dt = k2[AC] eq.26
d[AC]/dt = -k2K[C][AC]/(1 + K ([A] + [C])) eq.27
d[C]/dt = -d[AC]/dt eq.28
が得られる。これらの関係式は定常状態近似の成立しないような条件下でも適用できるので、Michelis-Mentenn式より広範なデータの解析に用いることが出来る。
なお、このような反応ではＡ、ＣおよびＡＣの初期濃度はＢの初期濃度を０として、 eq.19、eq.21、eq.22 から算出されるが、spana_redap
ではfunction2() としてこの初期濃度算出ルーチンを定義する。プログラム実行時、c[1], c[2] として[A]T, [C]Tを与え、[A], [C], [AC] の初期濃度として"99999"を設定すると、spana_redap はこのルーチンを使って各初期濃度を算出する。
A : x[1]、B : x[2]、C : x[3]、AC : x[4]、K : p[1]、k2 : p[2] 、[A]T : c[1]、[C]T : c[2] として、function関数を以下に示す。
/* Michaelis-Menten Reaction FOR REDAP */
/* A + C <-> AC -> B + C */
/* K = [AC]/[A][C] */
/* d[B]/dt = k[AC] */
/* p[1] = K , p[2] = k */
/* x[1] = [A] , x[2] = [B] , x[3] = [C] , x[4] = [AC] */
/* c[1]= [A]o + [AC]o , c[2] = [C]o + [AC]o */
#include <math.h>
extern double dx[],x[],c[],p[];
void function() /* Definition of kinetic processes */
{
dx[1]=-p[2]*x[4]*(1+p[1]*x[1])/(1+p[1]*(x[1]+x[3]));
dx[2]=p[2]*x[4];
dx[3]=p[1]*p[2]*x[3]*x[4]/(1+p[1]*(x[1]+x[3]));
dx[4]=-dx[3];
}
void function2() /* Initial concentration for equilibrium process */
{
double tem;
tem=p[1]*(c[1]+c[2])+1;
x[4]=(tem-sqrt(tem*tem-4*p[1]*p[1]*c[1]*c[2]))/2/p[1];
x[1]=c[1]-x[4];
x[3]=c[2]-x[4];
}
これらのプログラムでは解析する速度過程の対象は化学種濃度であるが、測定対象が電子スペクトルなど間接データである場合はスペクトルの時間経過式をさらに定義する必要がある。
この為には例えば、A、ＡＣ、B のモル吸光係数をそれぞれ εA、εAC、εB として吸光度 Abs の時間変化、
dAbs/dt = εAd[A]/dt + εACd[AC]/dt + εBd[B}/dt
を
dx[5]=p[3]*dx[1]+p[4]*dx[4]+p[5]*dx[2];
としてfunction内に追加定義する必要がある。
注
１） プログラムは通常Damping Factorの初期値を１とし、以降二乗残差の減少に伴って自動的に縮小・拡大しながら最適解の探索を行うが、SPANAでは系の非線形性が大きい場合など最小二乗計算画面から初期値を変更して探索を行うことも出来る。
２）本ＨＰでは、より一般的な非線形最小二乗法のプログラム、LSPE および REDAP、を公開している。"Download"のページを参照。
3) spana_lsq.c、および spana_redap.c は Microsoft Visual Studio C++ のC言語仕様で記述されているので、他のコンパイラ系では一部各仕様に合わせたソース・プログラムの書き換えが必要な場合がある。また、コンパイラはユーザー使用のWINDOWSが32bitか64bitかに依存して、それぞれシステム対応のコンパイラを使用しなくてはならない。
4）この速度解析用のプログラミングについては、以下に詳しく説明されている。
「生体機能関連化学実験法」、日本化学会生体機能関連化学部会編、p49、化学同人、 2003
