Go to the first, previous, next, last section, table of contents.


曲線近似の詳細

一般線形モデルへの変換

Calc の最小二乗近似エンジンは多変量線形モデルしか扱えません。 より正確に言えば、 a f(x,y,z) + b g(x,y,z) + c h(x,y,z) 書式の任意のモデルが扱えます。 ただし a,b,c はパラメータ、x,y,z は独立変数です (もちろん、3つに限らずいくつでもかまいません)。

簡単な多変量または多項式近似においては、 モデルをこの形式に変換するのは簡単です。 例えば、モデルが a + b x + c x^2 なら、 f(x) = 1, g(x) = x, h(x) = x^2 です。

他のモデルでは、 Calc はいろいろな代数操作を使って、問題を次の書式に変換しようとします。

Y(x,y,z) = A(a,b,c) F(x,y,z) + B(a,b,c) G(x,y,z) + C(a,b,c) H(x,y,z)

ただし Y,A,B,C,F,G,H は任意の関数です。 Calc は全てのデータ点について Y, F, G, H を計算し、 標準の線形近似を行って A, B, C の値を見出し、 次に方程式解法エンジンを使って A,B,C の各項を a,b,c について解きます。

驚くほど多くのモデルが、この一般書式に変換できます。 ここでは2つの例を見て、その働きを理解しましょう。 2つの独立変数と 2つのパラメータを持つべき乗モデル y = a x^b は、 次のように書替えられます。

y = a x^b
y = a exp(b ln(x))
y = exp(ln(a) + b ln(x))
ln(y) = ln(a) + b ln(x)

この式は望みの書式に対し次のようにマッチします。 Y = ln(y), A = ln(a), F = 1, B = b, G = ln(x) Calc はこのように、yx の対数を計算し、 AB を直線近似し、そしてそれらを解いて a = exp(A)b = B を求めます。

もう一つの興味深い例は、分配法則に従い展開できる「二次式」モデルです。

y = a + b*(x - c)^2
y = a + b c^2 - 2 b c x + b x^2

これは次のようにマッチします。 Y = y, A = a + b c^2, F = 1, B = -2 b c, G = x (因数 -2B の替わりに G の中に変換することもできました)。 (訳注: そしてさらに C = b, H = x^2 とマッチする。)

ガウス関数モデルはひどく複雑に見えます。 しかしよく見ると、 エクスポネンシャルの内部は2次式モデルに似ていることが判ります。 エクスポネンシャルは開いて Y に放り込んでやれば良いのです。

線形化不可能なモデルの取扱い

ガウス関数にバックグラウンド(定数)が加わったような、 すなわち [Gaussian] + d のような式は 一般線形書式に書替えられないモデルの一例です。 もしこのようなモデルで近似しなければならない場合、最善策は次のとおりです。 まず、モデルを線形化可能とするために(邪魔な)パラメータ群を定数群に置換え、 それから手動でその定数群を調整しながら近似操作を繰返すのです。 近似の具合を比べるには、グラフに可視化したり、 I a F が返す近似良好性指数を調べたり、 目的に適した何か他の方法を使ったりします。 注意:モデルによっては、線形化するのに複数のやり方があります。 ガウス+dモデルの場合、バックグラウンド d を定数に設定するか、 標準偏差 b と平均値 c を定数に設定すると線形化できます。

パラメータを定数に置換えたモデルで近似を行うのは簡単で、 まずそのパラメータ変数に適当な値をストアし、そして変数プロンプトの局面で、 パラメータ・リストからそのパラメータ群を削除します。

最後の手段は、fit よりもむしろ、 汎用の minimize 関数を使うことです。 どちらの関数も、式( χ^2 和)中のパラメータを調整して式を最小化する点で、結局同じです。 a F コマンドは、χ二乗和に関する特殊知識を利用して、 遙かに効率的なアルゴリズムを使うことができますが、 a N コマンドは力ずくで同じ事をします。

折衷案としては、 線形化の障害となっている少数のパラメータを見分けて minimize を使い、 残りのパラメータ群は fit 呼出しにまかせることです。 ここで最小化されるのは、 χ^2 値であって、xfit 関数の 5番目の結果として値が返されます。

minimize(xfit(gaus(a,b,c,d,x), x, [a,b,c], data)_5, d, guess)

ここで gaus はバックグラウンド付きガウシアンモデルを表し、 data はデータ行列を表し、 guessminimize が必要とする d の初期推定値を表します。 この演算は(minimize が使われた場合の常として)、 天文学的ではないにしろ、異常に遅いです。

線形化変換時の注意

I a F [xfit] コマンドは、 線形化変換を介する場合に少し面倒になります。 その結果の 2番めの項目に示されるベクトルは、 オリジナル・パラメータ a, b, c に対応するものではなく、 「変換パラメータ」 A, B, C のベクトルです。 (訳注: 3番めの)共分散行列は、その変換パラメータ群に関するものです。 4番め(訳注: 原文には fifth とあるが誤記であろう)の項目は 変換式のベクトルです。 変換パラメータ群がオリジナル・パラメータ群と等しいとき、 これは空ベクトル `[ ]' になります。 つまり A = a, B = b, ... のような場合であって、 例えば多項式近似のように、モデルが初めからオリジナル・パラメータ表記のままで 線形であった場合に成り立ちます。 もし(訳注: モデルを線形形式に変形する必要性から)パラメータ群を変換した場合、 ここには各オリジナル・パラメータへの変換式のベクトルが入ります。 変換式中での変換パラメータの名前ですが、 A`fitdummy(1)', B`fitdummy(2)', ... として表現されます。

Calc は結果を返すとき、 変換式中の `fitdummy(1)' に最初の変換パラメータの値を代入し、 他の変換パラメータの値も同様に代入して変換式群の値を求め、 これをオリジナル・パラメータ群の値としてオリジナルモデルに代入します。 H a FI a F の場合には、 オリジナル・パラメータの結果は誤差型式で算出する必要があるので、 Calc は共分散行列の対角要素群の平方根を変換パラメータ群の誤差値として採用し、 そこから Calc の標準誤差型式演算を開始します。

もしあなたが I a F で非線形モデルを使うなら、 共分散行列がオリジナル・パラメーターではなく、 変換パラメータに対応していることをよく認識しておいてください。 その共分散値を、変換式の存在下でどう解釈するか、それはあなた次第です。

入力が誤差型式を含む場合もちょっと込み入っています。 独立変数および従属変数 x, y, z があって、 そのうちひとつ以上が誤差型式でデータ入力されたとします。 Calc は誤差の二乗和のルートを取ることによって、 すべての誤差値を結合します。 そして xy を誤差なし数値に変えて、 z を結合済み誤差つきの誤差型式に変えます。 線形変換されたモデルの Y(x,y,z) 部分が計算され、 その結果は誤差型式となります。 その誤差部分は対応するデータ点の sigma_i に使われます。 もし何らかの理由で Y(x,y,z) が誤差型式の結果を返さなければ、 結合済み誤差(訳注: つまり z の誤差)が直接 sigma_i に使われます。 最後に、 F(x,y,z)G(x,y,z)などを計算する際には z の誤差も剥ぎ取られ、 線形変換されたモデルの右辺は普通の誤差無し演算で計算されます。

(これらの仕様は複雑に見えますが、 Y(x,y,z)が従属変数 zにのみ依存したり 単にzに等しいような、 実際によくある典型的なケースにおいて、 もっとも妥当な処理をするように設計されたものです。 多項式や多変量線形モデルのような普通のモデルでは、 結合済み誤差は単純にデータ点の sigma として使われ、それ以上の面倒はありません。)

線形化変換機構

使いたいモデルが線形化変換可能であるのに Calc の組込み規則がそれを理解できないこともあるでしょう。 Calc は代数書替え機構を使って線形変換します。 書替え規則は変数 FitRules の中に保存されていて、 s e FitRules コマンドを使うと、この規則を編集できます。 実は、FitRules を専用に編集する特別コマンド s F が存在します。 その他の変数操作 参照 。

書替え規則の議論は 書替え規則(Rewrite Rules) 参照 。

Calc は FitRules を次のように使用します。 まず、オリジナル・モデルを式形式に変換し、 モデル式を関数呼び出しfitmodel の引数とします (これは実際に定義された関数ではなく、書替え規則の目印になるだけです)。 パラメータ変数群は 関数呼び出し `fitparam(1)', `fitparam(2)', ... に リネームされ、独立変数群は `fitvar(1)', `fitvar(2)', ... に リネームされます。 従属変数は一番最後の fitvar に割りあてられます。 例えば、べき乗モデル a x^by = a x^b に変換され、 そして次のようになります。

fitmodel(fitvar(2) = fitparam(1) fitvar(1)^fitparam(2))

そして Calc はあたかも `C-u 0 a r FitRules' のように、 書替えを適用します。 (接頭引数のゼロは、それ以上進めなくなるまで書替え続けろ、という意味です。)

書替えが完了したとき、fitmodel 関数は次のように fitsystem 関数に置き換えられているはずです。

fitsystem(Y, FGH, abc)

ここで Y は関数 Y(x,y,z) を記述した式であり、 FGH は式 [F(x,y,z), G(x,y,z), H(x,y,z)] のベクトル、 abc はパラメータ変換式のベクトルで、 式中では変換パラメータ A`fitdummy(1)' として参照し B`fitdummy(2)', ... というように 参照します。 変換パラメータの数(FGH ベクトルの長さ)は オリジナル・パラメータの数(abc ベクトルの長さ)に 普通は一致しますが、これは必要条件ではありません。

べき乗モデルはやがて、次のように煮詰まります。

fitsystem(ln(fitvar(2)),
          [1, ln(fitvar(1))],
          [exp(fitdummy(1)), fitdummy(2)])

実際の FitRules の実行は複雑で、4段階で進行します。 まず普通の書替えを実施して、線形項を集めたり、 expln のような関数を分離して (Y に格納できるように)ずっと「外」に持っていったり、 あるいは(abcFGH に格納できるように)ずっと「内」に 集めたりします。 特に、定数でないべき乗は log/exp 形式に変換され、 カッコの掛け算は分配法則を使って展開されます。 割り算は `fitinv' 関数を使うように書き替えられます。 FitRules が稼働している間、`fitinv(x)'1/x を意味します。 (fitinv を使うと、線形的な項の認識が容易になります。) もし FitRules を改造するなら、 たぶんこの段階の規則群を修正するだけで済むと思います。

第2段階の規則は、実は第1, 第3段階でも適用されていて、 まず fitmodel を2引数形式の `fitmodel(Y, model)' に 書替えます。 ここで Y は 初期値ゼロであり、 modela=b から a-b 形式に変わっています (訳注: Y は左辺で model は右辺のように扱われる)。 次に model の表層にある逆変換可能な関数を引き剥がして それを Y に放り込もうとします。 逆関数の作成には方程式解法エンジンが呼び出されます。 最後に、逆変換がそれ以上不可能になったら、 fitmodel は 4引数の fitsystem に変えられます。 ここで、4番目の引数は model で、FGHabc の ベクトルは初期値が空です。 (最後のベクトルは今のところ、変換パラメータに対応して、本当に ABC です。)

第3段階は model の和の各項を `a*b*c' の意味する関数 `fitpart(a, b, c)' に変換します。 ここで、a は変数を含まない全ての因子、 b はパラメータだけを含む因子、 c は独立変数だけを含む因子です。 (もしこの分解が不可能ならば変換が完了せず、 Calc は "the model is too complex" と文句を言うでしょう。) 次に必要な変換パラメータの種類を最小限にするために、 bc 成分に等しい fitpart群が 分配法則を使ってまとめられます。

第4段階では、fitpart 項が FGHABC ベクトルの中に 移動されます。 さらに、式が効率的に計算できるように、 第1段階で実施された代数展開が今度は復旧されます。 最後に、もういちど方程式解法エンジンを呼び出して ABC ベクトルを abc ベクトルに変換します。 そして4番目の(この時点ではゼロになっている) model 引数が取り除かれ、 3引数の fitsystem が得られます。 これが最小二乗法エンジンの求める解です。

FitRules 関連で便利なコマンドに `hasfitparams(x)'`hasfitvars(x)' があり、 x がパラメータや独立変数を参照するかどうか調べます。 特に、もし引数がfitparam (あるいは fitvar)を 含んでいれば"true"を返し、そうでなければ"false"を返します。 ("true"は nonzero、"false"は zero にあたります。 実際に返される nonzero値は、式中に現れる全ての `fitparam(n)' あるいは `fitvar(n)' のうちの 最大の n が返されます。)

代数表記におけるデフォルト引数

代数表記における fit 関数は通常4つの引数を取ります。 `fit(model, vars, params, data)' ここで、modela F ' に続いて入力されるモデル式、 vars は独立変数またはそのベクトル、 params は同じくパラメータ群、 data はデータ行列です。 model を等式で与える場合、vars の長さ(変数の数)は data の行数と一致しなければならないことに注意してください。 model が等式を含まないベタ式なら 行数マイナス1 になります。 (実のところ、vars には従属変数の名前も許されますが、 ベタ式の場合は無視されます。)

params が省略された場合、model に現れる変数のうち、 vars を除いた全てがパラメータとなります。 vars までも省略されたら、 Calc はmodel に現れる全ての変数をアルファベット順にソートして、 前の方を params とし、後ろの方を vars とします。

別の方法では、長さ2〜3のベクトル modelvec によって 上記のモデルと変数を記述した、 `fit(modelvec, data)' が許されます。

もし近似ができなかった場合、fit 関数はシンボリック形式のまま 残されます。たいていの場合は注釈が付いて、 線形化変換に失敗した場合は"Model expression is too complex"と出ます。

efit (H a F に対応)や xfit (I a F)でも全く同様です。


Go to the first, previous, next, last section, table of contents.     利用度数