English version is here.
このブログで使われているプログラムのダウンロード・更新についてのガイドはこちら。
データを眺めた後でモデルを構築してそれがデータを説明していることを確かめる、というのはデータ解析の一つの理想的な終わり方のような気がします。しかし実際にやってみると作業が多くてなかなか骨が折れます。やる事はだいたい決まっているのに1つのツールではなかなかできず、プログラムも冗長になりがちです。今回の記事では非線形関数フィッティングをした上で更にその適合度(goodness-of-fit)を簡単に計算するツールを紹介します。
作ったモデルがデータに適合しているか確認するには、まず数式を用いてモデルを表現し、(多くの場合非線形の)その式をデータにフィッティングし、適合度の検定をします。以前 Curve Fitting Toolbox を持っていなかったときにこれらの作業をするのはとても困難でした。素のMatlabで非線形関数フィッティングに使える関数は ‘fminsearch’ だけでした、’fminsearch’はフィットしたパラメータの誤差の解析を行わないので、パラメータの解釈もそのパラメータに対する検定等も面倒になります。
なんとか非線形関数フィットが便利にできる関数がないものか探したところ、Fminuitというツールを見つけました。FminuitはCERNで開発されたMinuitをMatlab用にラッピングしたものです。詳しい説明は本家のページに譲りますが、かなり強力なツールです。このプログラムは数十年にわたってCERNでのデータ解析に活発に使われており、毎年ペタバイト単位で集まってくるデータを捌くのに使われています。Fortranを基礎とするコードは堅牢かつ高効率かつ多機能で、それでいてシンプルなコマンドで柔軟に使うことができます。本家でWindowsとLinux用のバイナリを配布しているので、私のGithubではmacOS用(Intel Mac)のみを置いておきました。Mac以外を使っている方はバイナリを本家からダウンロードしてくるか、コードを落としてきてコンパイルし、Matlabのパスに入れてください。
もちろん Fminuit の時点で既に強力なツールなのですが、以前紹介した ErrorNum クラスを用いて更に使い易くしようと試みました。’MinuitFitEN’ という関数を書いたのですが、実例を見る前に引数の確認のためヘルプを見てみましょう。
function res = MinuitFitEN(func, initp, x, y, varargin)
func - fit funciton handle func(p, x) % 関数ハンドルを引数として渡します。
initp - (nPar vec) estimated inital parameters % パラメータの初期値。
x - (any) 2nd argument of func % func の2番目の引数に該当するものをここに入れます。
y - (any) ErrorNum data to fit including error % フィット対象のデータを ErrorNum にして入れます。
-- from here optional parameters (default) % オプション群です。括弧内はデフォルト値です。
'Extradim' - (scalar) extra degree of freedom to subtract (0) % もし束縛条件などがあって、余計に引いておきたい自由度がある場合使います。
'LowerLimit' - (nPar vec) lower bound of the parameters (false) % パラメータの下限値がある場合に使います。
'UpperLimit' - (nPar vec) upper bound of the parameters (false) % 上限値も同様です。下限を設定したら上限も設定してください。
'StepSize' - (nPar vec) step sizes of the parameters (false) % フィッティングのステップサイズを指定したい場合に使います。
'minfunc' - (String) function to minimize ('PoissonNLL') % 最小化する関数です。デフォルトはポアソン分布です。カイ二乗フィットをするときは'Chi2'と書いてください。
'command' - (String) command to give to fminuit % fminuitに渡すコマンドを指定します。minuitのマニュアルを参照してください。
('set pri -1; scan; minimize; simplex; improve; minos')この関数を使うと非線形関数フィットが適合度検定も含めて比較的簡単にできるようになります。まずはフィットするデータを準備します。
x = 0:0.2:3;
n_obs = 30; % 測定の回数に該当するものです
sigma = 5; % データの標準偏差
% フィッティングを用いて以下のパラメータを推定します
base = 2; % 指数関数の底
offset = 3; % y軸のオフセット
% データを構築します
% 最後の項で測定誤差を加えています。
y_obs = repmat(base.^x, n_obs, 1) + offset + sigma * randn(n_obs, length(x));
y_err = ErrorNum.create(y_obs, 1); % ErrorNum.create は 標準誤差を計算します。
figure(101);
plotEN(x, y_err, '.', 'capsize', 0); % 最後に EN とつく関数は ErrorNum 用の関数です。
これを実行すると、このようなデータが生成されます。だいたい指数関数ですが、測定誤差(今回のは人為的に入れられたものですが)のせいで値が上下にふらついています。これに対して関数フィッティングをやってみましょう。
powfunc = @(p, x) p(1).^x + p(2); % これがモデル関数です。
% パラメータに適当な初期値を与えます。真値に近いほどフィッティングが成功しやすくなります。
base_guess = 3;
offset_guess = 1;
initial_p = [base_guess, offset_guess];
fitresult = MinuitFitEN(powfunc, initial_p, x, y_err, 'minfunc', 'Chi2');最初の行で無名関数を作ってハンドルを手に入れています。もしモデル関数用の.mファイルが既にある場合第一引数を ‘@関数名’ としてハンドルを渡してください。フィッティングに必要なコードはこれだけです。 MinuitFitEN の最後の引数2つはフィッティングのタイプ(最小化する量)を示しています。今回は誤差がガウス分布なのを想定しているので、カイ二乗フィット’Chi2’を選択しています。私が個人的にニューロンの発火数に対する解析をよくやるので、デフォルトの設定がポアソン尤度関数’PoissonNLL’になっていますが、ANOVA (分散分析)が適用可能な解析の場合は’Chi2’の方を選択してください。(‘Chi2’の方がデフォルトとして適切な気もするのでこちらをデフォルトに変更するかもしれません。)
さて、フィッティングの結果を見てみましょう。
>> fitresult.p
ans =
ErrorNum:
2.07 ± 0.06, 2.8 ± 0.3MinuitFitEN は、推定されたパラメータ(p)を含む構造体を返します(誤差もちゃんと計算されています)。この誤差はANOVAにおける 1 標準偏差なので、だいたい68%の確率で真値がこの誤差の範囲に入っているはずです。さきほど設定した底の値は 2 なので、少しだけ外れています。オフセットは 3 なので誤差の範囲に入っています。どうやらちゃんと動いているようです。パラメータ推定をするときにいつも誤差を計算してくれるのは安心感があります。出力が ErrorNum なので、桁数の調整なども必要なく、適切な有効数字で返してくれます。
さて、先ほども述べた通り MinuitFitEN は関数の適合度も計算してくれます。
>> fitresult.redchi2
ans =
1.5140
>> fitresult.redchi2sig
ans =
0.9033‘redchi2’ は reduced chi-squared の略で、カイ二乗を自由度で割った値になります。詳しい説明はWikipediaに任せますが(英語版になります。日本語の説明はありませんでした…。)、この値が1に近ければフィッティングがうまくいっているということになります。小さすぎるときはオーバーフィット(あるいは誤差の過大評価)、逆に大きすぎるときはアンダーフィット(あるいは誤差の過小評価)していると解釈できます。
1に近いといっても、どのくらい近ければ許されるのか明確ではないので、このカイ二乗の値をカイ二乗分布と比較して、どの程度稀なオーバーフィットやアンダーフィットなのかを調べられるようにしたのが ‘redchi2sig’ の値になります。これは適合度を評価する際のp値のようなものです。仮にこの値が0.99の場合は、フィッティングがうまくいっていると仮定して1%程度の確率でしか得られない値になっている(したがって、アンダーフィット(モデリングがうまくいっていない)の可能性が高い)、というように解釈します。モデリングがうまくいっていることを確かめるためにはこの値があまり極端な値をとらないことを示すことが必要になります(状況にもよりますが、0.05 から 0.95 あたりの間に入っているのが目安かと思います。p値を使うときと同じような注意が必要になります。)
これらの情報はフィッティングがうまくいっていることを確認するために大変有用な情報で、私自身も過去に書いた論文(Ito et al., 2017; Ito et al., 2020)で積極的に使用しています。しかし、敷居の高さもあってか一部の学問領域を除いて浸透しているとは言い難い状況です。これから書く予定の記事の中でも、これらの統計は使用していくことになると思いますので、そのときにまた具体例から適切な活用方法を掴んでいただけると幸いです。例に漏れず、これらの関数も隠れた高度な使い方がたくさんありますので、それらも紹介していけると嬉しいです。もし使ってくださる奇特な方がいらっしゃる場合は、それまでヘルプを参照して使用方法を確認していただけると有り難いです。
