odatse_plt_model_evidence#

機能概要#

model evidence を計算します

書式#

odatse_plt_model_evidence [OPTION]... -n NDATA FILE...

説明#

PAMC の出力ファイル FILE から beta および分配関数の値を取り出し、model evidence を計算します。計算結果は標準出力およびテキストファイル(デフォルト: model_evidence.txt)に出力されます。また、結果をプロットした図を画像ファイル(デフォルト: model_evidence.png)に出力します。

複数の FILE を指定した場合は、それらの model evidence の平均と標準偏差を出力し、エラーバー付きのプロットを生成します。

注釈

ODAT-SE 本体と同様に Python 3.9 以上が必要です。
計算はすべて対数スケールで行われるため、数値的に安定しています。
プロットのx軸(beta)は常に対数スケールで表示されます。

FILE: PAMCの出力ファイル名 (fx.txt)。複数のファイルを指定可能です。
-n NDATA, --ndata NDATA: 各データセットに含まれるデータ点の数をカンマ区切りの整数値で指定します。必須パラメータです。例: 「100」(1つのデータセットに100点)、「50,100,75」(3つのデータセットにそれぞれ50点、100点、75点)
-w WEIGHT, --weight WEIGHT: データセット間の相対重みをカンマ区切りの数値で指定します。重みは和が 1.0 になるように自動で規格化されます。重みの数値の個数とデータ点の個数は一致させる必要があります。
-V VOLUME, --Volume VOLUME: 事前確率分布の規格化因子 (定義域の体積 $V_\Omega$) を指定します。デフォルトは 1.0 です。
-f RESULT, --result RESULT: model evidence の値を出力するファイル名を指定します。デフォルトは model_evidence.txt です。
-o OUTPUT, --output OUTPUT: model evidence のプロットを出力するファイル名を指定します。出力形式は拡張子を元に設定され、matplotlib がサポートする形式を指定可能です。デフォルトは model_evidence.png です。
--auto-focus: model evidence の最大値をもとに、プロットの表示範囲を自動決定するフラグです。範囲の絞り込みの強さは --focus-factor オプションで制御します。
--focus-factor: --auto-focus の絞り込みの強さを指定します。0 から 1 の実数（小さいほど狭い）で、--auto-focus が指定されていない場合は無効です。デフォルトは 0.5 です。
-h, --help: ヘルプメッセージを表示してプログラムを終了します。

使用例#

基本的な使用方法(1つのデータファイルと1つのデータセット)
```
$ odatse_plt_model_evidence -n 100 fx.txt
```
データ点数100のデータセットについて model evidence を計算し、 model_evidence.txt と model_evidence.png を出力します。
複数のデータセットがある場合
```
$ odatse_plt_model_evidence -n 50,100,75 -w 0.2,0.5,0.3 fx.txt
```
3つのデータセット(データ点数がそれぞれ50、100、75で、相対重みが0.2、0.5、0.3)について、 model evidence を計算します。
複数のデータファイルを使用する場合
```
$ odatse_plt_model_evidence -n 100 -o evidence_plot.pdf -f evidence_data.txt fx_1.txt fx_2.txt fx_3.txt
```
3つのデータファイルから model evidence を計算し、それらの平均と標準偏差を求めます。結果を evidence_data.txt に出力し、evidence_plot.pdf にエラーバー付きのプロットを生成します。

補足事項#

model evidence の計算#

R-factor を以下のように定義します:

\[R(X;D)^2 = \sum_\mu w_\mu \sum_i \left( I_\mu(\theta_i) - I^{\text{(cal)}}_\mu(\theta_i;X) \right)^2\]

$I_\mu(\theta_i)$ はデータセット $\mu$ に含まれる計測データ点、 $I^{\text{(cal)}}_\mu(\theta_i;X)$ はパラメータ $X$ のもとでの理論計算値とします。 $w_\mu$ はデータセットの相対重みで、合計が1になるよう規格化されています。

model evidence $P(D|\beta)$ は以下の式で計算されます:

\[\log P(D|\beta) = \log Z(D;\beta) - \log V_\Omega + \frac{n}{2} \log \beta + \sum_{\mu} \frac{n_{\mu}}{2} \log w_{\mu} - \frac{n}{2} \log \pi\]

ここで、 $Z(D;\beta)$ は分配関数

\[Z(D;\beta) = \int_\Omega \exp\left(-\beta\,R(X;D)^2\right) dX\]

また、

$V_\Omega$: 事前確率分布の規格化因子
$n_\mu$: 各データセットのデータ点数
$n$: 全データ点数(すべてのデータセットの合計)
$\beta$: 逆温度

とします。

入力ファイルの形式#

入力ファイル(PAMCの出力ファイル)は以下の形式を想定しています:

# コメント行(任意)
beta_value  fx_mean  fx_var  nreplica  logz_value  acceptance
...

スクリプトは各行から以下の値を読み取ります:

第1列(インデックス0): beta値(逆温度)
第5列(インデックス4): logz値(対数分配関数)

出力ファイルの形式#

出力ファイル(model_evidence.txt)の形式は以下の通りです:

# max log_P(D;beta) = {最大値} at Tstep = {インデックス}, beta = {対応するbeta値}
# $1: Tstep
# $2: beta
# $3: model_evidence
0  beta0  model_evidence0
1  beta1  model_evidence1
...

複数の入力ファイルを処理した場合は、標準偏差の列が追加されます:

# max log_P(D;beta) = {最大値} at Tstep = {インデックス}, beta = {対応するbeta値}
# $1: Tstep
# $2: beta
# $3: average model_evidence
# $4: standard deviation
0  beta0  avg_model_evidence0  std0
1  beta1  avg_model_evidence1  std1
...

処理の仕組み#

このスクリプトは以下の手順で処理を行います:

入力ファイルからbeta値とlogz値を読み込みます
各データセットのデータ点数と重みを取得します
model evidenceの対数値を計算します
複数ファイルの場合は平均と標準偏差を計算します
結果をファイルに出力します
model evidenceをbetaの関数としてプロットします

プロットの特性#

X軸(beta)は常に対数スケールで表示されます
単一ファイルの場合は点のみ、複数ファイルの場合はエラーバー付きで表示されます
マーカーは赤色の「x」で表示されます
グリッド線を表示し、データの位置を把握しやすくします

エラー処理#

入力ファイルが存在しない場合: ファイルオープンエラーが発生します
データ形式が不正: numpy.loadtxtでエラーが発生します
NDATAとWEIGHTの長さが一致しない場合、またはデータ点数・重みが正でない場合: エラーメッセージを表示して終了します

特に、データ点数のリストとそれらの重みの数は必ず一致させる必要があります。