Pythonでスペクトルのフィッティング

非プログラマーがPythonを学ぶために読んだ本 - 最終防衛ライン3 の「2. 一ヶ月勉強して、できるようになったこと」で書いた、Pythonでスペクトルのフィッティングを行った例を示します。matplotlib と pandas などを扱ったグラフの出力に関しては、別当まとめようと思います。matplotlib は手に取ってみるとなんとなくで扱えるのですが、真面目に取り組むと結構ややこしく、自分のためにも備忘録をしっかりまとめておきたい。
フィッティングだけならExcelでもできますが、データを大量に取り扱うとなるとマクロを組む羽目になりますし、マクロよりも、スクリプト言語を扱った方が拡張性があります。デーだ処理だけなら、Rを使う手もありますが、グラフを出力すること以外にもPythonを使った方が、やりたいこと、またできることが多いので。

フィッティングする適当なデータがなかったので、放射線スペクトル表示ツール SPViewer / 説明ページ における スペクトル例(FNF-401) を利用しました。確認用データを生成してそれをフィッティングしても良いのですが、それだとフィッティングできて当たり前ですし。
放射線スペクトルは、正規分布を取ることが知られているのでガウシアンフィッティングに最適。またバックグラウンドも二次式程度で近似できそうなので、丁度良いかなと。

f:id:lastline:20180412092647p:plain

スペクトルのフィット:例2Pythonでやった結果を表示しています。
先ず、pandas で CSV ファイルを読み込み、
フィッティング関数を定義し、scipy.optimize の curve_fit でフィッティングを実行します。
そして、その結果を matplotlib で描画しています。
フィッティング関数の定義や、データの取り扱いには numpy を利用しています。

以下のコードを実行すると、グラフが出力されます。
CSVファイルは、先述のように 放射線スペクトル表示ツール SPViewer / 説明ページ から スペクトル例(FNF-401) をダウンロードし、同名のフォルダに保存してください。URLから直接取得することもできますが、本質ではないのでカット。
過剰とも思えるコメントを付記しています。調べた際に、分かりにくかったポイントをまとめてあります。。そのくせに、関数の document は記述していない。

#モジュールのインポート
#import モジュール名 as * は非推奨。使っているコードを見かけますけど。
import pandas as pd
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

#データの読み込み
#関連として df を用いるが、データに即した変数名を使うべき。
#ベクトルとして取り出すことも可能。
#ただし、x, y を個別に取り出した方がフィッティングの際に扱いやすい。
#engine='python' と指定すると日本語ファイルも読み込めます。
file = 'fnf401_milk.csv'
df = pd.read_csv(file, engine='python', skiprows=250, nrows=350, header=None)

#元データがint値で、そのままではフィッティングできないので、float型に変換している。
channels = np.array(df.iloc[:, 0], dtype=float)
counts = np.array(df.iloc[:, 1], dtype=float)

#フィッティング関数の定義
#n次式を定義しやすいように、0次から記述した。
def quadratic(x, a, b, c):
    return a + b * x + c * x**2

#ガウシアンはより簡素に記述できるが、初期値を類推しやすい式にした。
def gauss(x, a, sigma, mean):
    return a / np.sqrt(2.0*np.pi) / sigma * np.exp(-((x-mean)/sigma)**2/2)

#フィッティング関数は二次関数とガウシアンをまとめて定義してもよい。
#ガウシアンを三度呼び出すので、分けて記述した。
#文頭のアスタリスクのついた*paramaters は引数を幾らでも渡せる。
#本来ならば引数を個別に指定すべき。
#今回はあまりにも冗長になるため、*paramaters と指定した。
def fitting(x, *paramaters):
    a0, b0, c0, a1, sigma1, mean1, a2, sigma2, mean2, a3, sigma3, mean3, = paramaters
    return quadratic(x, a0, b0, c0) + \
           gauss(x, a1, sigma1, mean1) + \
           gauss(x, a2, sigma2, mean2) + \
           gauss(x, a3, sigma3, mean3)

#初期値の設定。
#ガウシアンの場合、meanを適切に設定しないと発散しやすい。
#線形近似であっても適切に設定することが望ましい。
#関数の引数が複数ある場合は、タプルで渡すこと。
initial_quadratic = 70., 0.1, 0.1 #a0, b0, c0
initial_gauss1 = 1., 1., 300. #a1, sigma1, mean1
initial_gauss2 = 1., 1., 325. #a2, sigma2, mean2
initial_gauss3 = 1., 1., 395. #a3, sigma3, mean3

#初期値が多いので、個別に指定して接続した。
initial_parameter = initial_quadratic + initial_gauss1 + initial_gauss2 + initial_gauss3

#フィッティングの実行
#paramater_optimalが最適化されたパラメータで、covarianceは共分散。
paramater_optimal, covariance = curve_fit(fitting, channels, counts, initial_parameter)

#matplotlibの設定
plt.rc('font', family='Arial', size=14) #フォントの設定
plt.rc('xtick.major', width=1, size=6) #x軸の主目盛りの設定
plt.rc('xtick', direction='in', top=True) #x軸目盛りの向き、上側に表示するか
plt.rc('ytick.major', width=1, size=6) #y軸の主目盛りの設定
plt.rc('ytick', direction='in', right=True) #y軸目盛りの向き、右側に表示するか
plt.rc('axes', linewidth=1.5) #枠線の太さ
plt.rc('lines', linewidth=1.0) #プロットの太さ

#今回は、Axes が一つしかないので fig, ax1 = plt.subplots() の方が簡素。
#慣例として ax1 を用いるが、df 同様に意味のある変数名にすべき。
fig = plt.figure(figsize=(8, 6), dpi=80)
gs = gridspec.GridSpec(1, 1)
ax1 = fig.add_subplot(gs[0])

#データをグラフに表示する。
ax1.plot(channels, counts, ls='-', c='g', label='row data')
ax1.set_xlabel('Channel', fontsize=16)
ax1.set_ylabel('Count', fontsize=16)
ax1.yaxis.set_data_interval(0, 100) #set_ylim を使う方が一般的

#フィッティングした曲線の表示
#np.corrcoef で相関係数を求める。
#結果が2×2の行列で返るので、1行目2列目もしくは2行目1列目の値を取得。
ax1.plot(channels, fitting(channels, *paramater_optimal),
         lw=1.5, ls='--', c='k', label='fitting curve\n' + \
         f'$R^2$ = {np.corrcoef(counts, fitting(channels, *paramater_optimal))[0][1]}')

#二次関数のパラメーターを利用したいので、スライスで取り出す。
a0, b0, c0 = paramater_optimal[0:3]

ax1.plot(channels, quadratic(channels, a0, b0, c0), c='b',
         label=r'$y = a_0 + b_0x + c_0x^2$' +'\n'+ f'$a_0$ = {a0}\nb$_0$ = {b0}\nc$_0$ = {c0}')

#各meanをグラフ上に表示する。繰り返しのため for 文を利用。
for i in range(3,12,3):
#この関数内の、a0, b0, c0 はグローバル変数である。
    def quadratic_gauss(x, *paramaters):
        a1, sigma1, mean1, = paramaters
        return quadratic(x, a0, b0, c0) + \
               gauss(x, a1, sigma1, mean1)

#annotate で Axes 内にテキストや矢印を表示できる。
    mean = paramater_optimal[i+2]
    ax1.annotate(int(mean), fontsize = 14,
                 xy=(mean, quadratic_gauss(mean, *paramater_optimal[i:i+3])+3),
                 xytext=(mean, quadratic_gauss(mean, *paramater_optimal[i:i+3])+15),
                 ha="center", va='top',
                 arrowprops=dict(shrink=0.05, facecolor='k'))

#凡例の表示。Figure の判例を出力することもできる。
ax1.legend(loc='upper right', frameon=False)

gs.tight_layout(fig) #グラフ描画前に行うべし。慣例では、fig.tight_layout()
plt.show(fig) #慣例では plt.show() と fig を指定しない。

参考文献やサイト

詳細! Python 3 入門ノート

詳細! Python 3 入門ノート

matplotlib や numpy が簡単ではありますが、まとまっています。

matplotlib や numpy の取り扱いで大変役に立ちました。

Python for Data Analysis: Data Wrangling with Pandas, NumPy, and IPython

Python for Data Analysis: Data Wrangling with Pandas, NumPy, and IPython

主に、pandas の取り扱いを参考にしました。