実験データの解析やレポートで、プログラムを少し真面目に作る

記事の内容は私個人の見解であり、所属する学科組織またはサークルを代表するものではありません.

はじめに

PhysiKyuアドカレ15日目です。九大理物B3の北川が担当します。Physikyuにも一応IT班的な下部組織あるのですが、一応その言い出しっぺである僕がここでIT系の話題を書いて、いい感じにPhysiKyuのフレッシュなメンバーにIT班を認知してもらいたいと思っています。

今回は、「データのプロットとフィッティングをするプログラム」という例を通して、

• 動くプログラムをどう組み立てていくか
• プログラムの再利用性

の2点に絞って伝えていきたいと思います。

良いプログラム

早速コードを書いていきたいところですが、その前に、コードを書くときに心がけておきたいことを話しておきます。僕の知っている限り、多くの物理学徒は、一般教養の授業でプログラミングを少しだけやっていて、長く使うプログラムに意識を向けたことがないと思います。ただ、ソフトウェア業界だと、出来るだけ多くの場面で支えて、多くの人が参加し、出来るだけ長く使えるプログラムを作っていきたいものです。細かいことを言えばキリがないですが、ともかくそういったプログラム作成のために良いコードには最低限次の3つの性質が備わっていると思います。

• 可読性
• 可変性
• 再利用性

もちろん、これらは完全に独立しているわけではないです。「可読性」は、あるコードがどういう処理をするものなのかがわかりやすいこと。「可変性」は、コードに変更を加えなければならないとき、変更箇所が出来るだけ少なく済むこと。「再利用性」は、コードがいろんな状況で使えることです。プログラマ界隈では、これらのために長い年月をかけていろんなノウハウが積み重ねられてきています。

動くプログラムを組み立てていく

それでは、コードを書いていきます。Pythonを使いますが、インストールとか環境構築とかについてはすっ飛ばします。

こんな状況を想定してみてください。

あなたは物理学科の学生です。必修実験がありますが、指定されたソフトがM1以上のMacでは動かしにくいことがわかりました。そこでMacユーザーのあなたは便利なソフトが使えず、PythonかExcelで解析するよう教員から言われます。Pythonでプログラムを書いておくと後々便利なので、勉強がてらPythonを使うことにしました。ただし、あなたはPythonでのデータ解析はほとんどやったことがないとします。( ほぼ実話です。)

今回ほしいプログラムは、「配布されたCSVファイルを読み込み、散布図にデータをプロットして、理論的に導出される式でフィッティングをする。また、フィッティング精度も計算する。」というものです。

まずはCSVファイルの読み込みから、とやりたくなるかもしれませんが、まずはこのプログラムの骨組みを見極めて、いったんゴールまで1本の道を通してしまうことが大切です。

1. データを用意する
2. データをフィッティングする
3. フィッティングの精度を評価する
4. データをプロットする

まずは1~3までをひとつなぎにできるように、本来やりたい処理を、実装が簡単な試験的な処理に置き換えてコードを書いていきます。今回の例だと、次のように置き換えます。

本来欲しい処理	簡単な処理
CSVを読み込む	y=sin(x)にノイズを乗せたデータを用意する
適切な関数でフィッティングする	y=Asin(Bx+C)でフィッティングする
決定係数で精度評価	残差2乗平均で評価
データをプロットし、タイトルや軸ラベルなどをつける	プロットするだけ

なんでこんなことをするのかというと、いきなり複雑な処理を書いていくと、実際にプログラムを実行できるようになるまで時間がかかり、動くかどうかわからないままコードを書いていかなければなりません。また、動かなかった場合、エラーの対応がやや煩雑になります。

以上の簡単な処理を、インターネットで調べながらPythonで書いてみます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy import stats

# y = sin(x)にノイズを乗せたデータを用意する
np.random.seed(0)
x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x) + np.random.normal(scale=0.2, size=x.shape)

# フィッティングする
def fitting_func(x, a, b, c):
    return a * np.sin(b * x + c)
params, params_covariance = curve_fit(fitting_func, x, y, p0=[1, 1, 0])

# 精度を評価する
r = y - fitting_func(x, *params)
print(f"残差2乗平均: {np.sum(r**2)/len(x):.3}")

# プロットする
plt.scatter(x, y)
plt.plot(x, fitting_func(x, *params))

plt.show()

結果はこんな感じになります。

残差2乗平均: 0.0348

悪くないですね。

これで一旦、大まかな流れが正しいことがわかったので、少し処理を複雑にしてみましょう。

本来欲しい処理	簡単な処理から、本来欲しい処理に少し近づけた処理
CSVを読み込む	y=sin(x)にノイズを乗せたデータが入ったCSVファイルを事前に作成して、pythonで読み込む
適切な関数でフィッティングする	y=Asin(Bx+C)でフィッティングする
決定係数で精度評価	残差2乗平均で評価
データをプロットし、タイトルや軸ラベルなどをつける	データをプロットし、タイトルや軸ラベルなどをつける

ちなみに、このCSVファイルの事前準備もPythonでサクッとできてしまうので、本筋とは離れますが一応共有です。

import numpy as np
import pandas as pd

# データ生成
np.random.seed(0)
x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x) + np.random.normal(scale=0.2, size=x.shape)

# データフレームの作成
data = pd.DataFrame({'x': x, 'y': y})

# CSVファイルに保存
csv_file_path = 'sin_data_with_noise.csv'
data.to_csv(csv_file_path, index=False)

話を戻します。先ほどの簡単な処理から、少しだけ本来欲しい処理に近づけて、また動作確認します。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit

# CSVファイルを読み込む
csv_file_path = 'sin_data_with_noise.csv'
data = pd.read_csv(csv_file_path)

# データを用意する
x = data['x']
y = data['y']

# フィッティングする
def fitting_func(x, a, b, c):
    return a * np.sin(b * x + c)
params, params_covariance = curve_fit(fitting_func, x, y, p0=[1, 1, 0])

# 精度を評価する
r = y - fitting_func(x, *params)
print(f"-"*30)
print(f"残差2乗平均: {np.sum(r**2)/len(x):.3}")
print(f"-"*30)

# プロットする
plt.scatter(x, y, label='observed data')
plt.plot(x, fitting_func(x, *params), label=f'y={params[0]:.3}sin({params[1]:.3}x+{params[2]:.3})', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Fitting y = sin(x) to observed data')
plt.legend()

plt.show()

結果はこんな感じになります。
------------------------------
残差2乗平均: 0.0348
------------------------------

いい感じになってきましたね。さらに目的の処理に近づけていきましょう。

本来欲しい処理	簡単な処理から、本来欲しい処理にさらに近づけた処理
CSVを読み込む	y=sin(x)にノイズを乗せたデータが入ったCSVファイルを事前に作成して、pythonで読み込む
適切な関数でフィッティングし、精度を評価する	y=Asin(Bx+C)でフィッティングし、精度を評価する
決定係数で精度評価	決定係数で精度評価
データをプロットし、タイトルや軸ラベルなどをつける	データをプロットし、タイトルや軸ラベルなどをつける

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit

# CSVファイルを読み込む
csv_file_path = 'sin_data_with_noise.csv'
data = pd.read_csv(csv_file_path)

# データを用意する
x = data['x']
y = data['y']

# フィッティングする
def fitting_func(x, a, b, c):
    return a * np.sin(b * x + c)
params, params_covariance = curve_fit(fitting_func, x, y, p0=[1, 1, 0])

# 精度を評価する
residuals =  y - fitting_func(x, *params)
rss = np.sum(residuals**2) 
tss = np.sum((y-np.mean(y))**2)
r_squared = 1 - (rss / tss)
print(f"-"*30)
print(f"決定係数: {r_squared:.3}")
print(f"-"*30)


# プロットする
plt.scatter(x, y, label='observed data')
plt.plot(x, fitting_func(x, *params), label=f'y={params[0]:.3}sin({params[1]:.3}x+{params[2]:.3})', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Fitting y = sin(x) to observed data')
plt.legend()

結果はこんな感じになります。
------------------------------
決定係数: 0.938
------------------------------

ここまできたら、あとは実際のデータを読み込み、適切なフィッティング関数に書き換えれば、やりたいことは実現できるわけです。この作業は実際に自然科学的なプロセスが必要なので、今回は割愛します。

再利用性に気を配る

以上で、今回の目的である「配布されたCSVファイルを読み込み、散布図にデータをプロットして、理論的に導出される式でフィッティングをする。また、フィッティング精度も計算する。」というプログラムは作成できました。が、ここで終わらないのが今回の記事のミソです。　

今回作成したプログラムを使い捨てにして、別のデータに対してまた同じようなプログラムを作るのはちょっと面倒です。
1.データを用意する, 2.データをフィッティングする, 3.フィッティングの精度を評価する, 4.データをプロットする, というプロセスは、今回配布されたCSVファイルだけに適用するものではありませんよね。

さらに、例えばフィッティング関数がよくわかっていない場合を想像すると、一旦フィッティングとかをすっ飛ばしてプロットしてみるという作業も必要ですね。

このように、一度作ったプログラムを、出来るだけ使いまわしたいし、いろんな場面で使えるように育てていきたいわけですが、ソフトウェア業界ではそのための概念と手法が、さまざまな試行錯誤の歴史を通してかなり発達しています。

この記事では、それらの手法を身につける前の、基本中の基本となる概念のうちの2つ、 関数と クラスを用いて、今回作成したプログラムをもう少し再利用性の高いものにしていきます。

関数

まずは関数ですが、関数については、コードの一連の流れにおいて、「サブタスクを実現しているコードをまとめておく」ものと思えば、とりあえず良いでしょう。ここでいっている「サブタスク」とは、今回の例で言うと、「CSVを読み込む」とか、「精度を評価する」とかいった、一連のコードの中でも論理的に仕事として切り分けて考えることができる部分のことです。

それでは、「CSVファイルを読み込んで、プロットできるデータを準備する」仕事をしている部分を関数にしてやりましょう。

def read_csv(file_name):
    # CSVファイルを読み込む
    csv_file_path = file_name
    data = pd.read_csv(csv_file_path)

    # データを用意する
    x = data['x']
    y = data['y']

    return x, y

このような関数をコードのどこかで定義しておいて、コードの一連の流れの中ではこの関数を使うようにします。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit

def read_csv(file_name):
    # CSVファイルを読み込む
    csv_file_path = file_name
    data = pd.read_csv(csv_file_path)

    # データを用意する
    x = data['x']
    y = data['y']

    return x, y

# ↓↓↓↓↓↓↓↓↓↓実際に実行される部分↓↓↓↓↓↓↓↓↓↓

# CSVデータを読み込み、データを用意
x, y = read_csv('sin_data_with_noise.csv')

# フィッティングする
def fitting_func(x, a, b, c):
    return a * np.sin(b * x + c)
params, params_covariance = curve_fit(fitting_func, x, y, p0=[1, 1, 0])

# 精度を評価する
residuals =  y - fitting_func(x, *params)
rss = np.sum(residuals**2) 
tss = np.sum((y-np.mean(y))**2)
r_squared = 1 - (rss / tss)
print(f"-"*30)
print(f"決定係数: {r_squared:.3}")
print(f"-"*30)


# プロットする
plt.scatter(x, y, label='observed data')
plt.plot(x, fitting_func(x, *params), label=f'y={params[0]:.3}sin({params[1]:.3}x+{params[2]:.3})', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Fitting y = sin(x) to observed data')
plt.legend()

同様にして、他のいくつかのサブタスクも関数にしてしまいましょう。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit

def read_csv(file_name):
    # CSVファイルを読み込む
    csv_file_path = file_name
    data = pd.read_csv(csv_file_path)

    # データを用意する
    x = data['x']
    y = data['y']

    return x, y

def calculate_fitting_accuracy(x, y, fitting_func):
    residuals =  y - fitting_func(x, *params)
    rss = np.sum(residuals**2) 
    tss = np.sum((y-np.mean(y))**2)
    r_squared = 1 - (rss / tss)

    return r_squared

def print_fitting_accuracy(indicator_name, value):
    print(f"-"*30)
    print(f"{indicator_name}: {value:.3}")
    print(f"-"*30)

def evaluate_fitting_accuracy(x, y, fitting_func):
    r_squared = calculate_fitting_accuracy(x, y, fitting_func)
    print_fitting_accuracy("決定係数", r_squared)


# ↓↓↓↓↓↓↓↓↓↓実際に実行される部分↓↓↓↓↓↓↓↓↓↓

# CSVデータを読み込み、データを用意
x, y = read_csv('sin_data_with_noise.csv')

# フィッティングする
def fitting_func(x, a, b, c):
    return a * np.sin(b * x + c)
params, params_covariance = curve_fit(fitting_func, x, y, p0=[1, 1, 0])

# 精度を評価する
evaluate_fitting_accuracy(x, y, fitting_func)

# プロットする
plt.scatter(x, y, label='observed data')
plt.plot(x, fitting_func(x, *params), label=f'y={params[0]:.3}sin({params[1]:.3}x+{params[2]:.3})', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Fitting y = sin(x) to observed data')
plt.legend()

このようにしてサブタスクを切り分けておけば、例えば、「フィッティングをすっ飛ばしていったん読み込んだデータをプロットしてみる」ようなシナリオも、コードを読んで必要なところをコピペしたり、調べ直したりするまでもなくすぐに対応できます。

x, y = read_csv('sin_data_with_noise.csv')
plt.scatter(x, y)
plt.show()

結果

「そんなのコピペすればいいじゃん」と思うかもしれませんが、それは後々自分を苦しめることになります。例えば今回の場合、「CSVデータを読み込み、データを用意する」というサブタスクを、関数にするのではなくコピペして対応したとしましょう。そして、この、「CSVデータを読み込み、データを用意する」というサブタスクを実行するロジックを少し変更したいとします。この場合、コピペしてしまったので、修正しなければならない箇所は２箇所になります。同じコードをコピペすればするほど、後々修正しなければならない箇所が増えていくわけです。関数として定義しておけば修正箇所は１箇所で済みます。

こんな感じで、関数を利用したい気持ちが伝わったでしょうか。

クラス

次に、これは中々イメージが湧きにくいし、正直今回のような短いコードではその威力がわかりづらいと思うのですが、クラスという概念を説明します。関数はいわば「機能を表す」ためのものでしたが、クラスはより抽象的に、「役割を表す」ものになります。「一連の処理の中で、何らかの一貫した役割を果たす"もの"を定義する」ものがクラスです。クラスは、取りうる状態と取りうる動作を規定してやることで定義できます。クラス自体は雛形みたいなもので、クラスに、状態を表す実際のデータを注入すると、クラスのインスタンスが生まれます。

今回作成したような短いプログラムでいろんなクラスを定義するのは難しいですが、プログラム自体が、「データを解析する」という役割を持ったものだと考えると、このクラスには、「解析すべきデータ」と、「フィッティング関数」、「フィッティング後のパラメータ」といった取りうる「状態」(ステート)があり、「データを読み込む」、「フィッティングをする」、「精度を評価する」、「プロットする」といった取りうる「動作」(メソッド)があります。このようなクラスをPythonで表現すると次のようになります。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit

class DataAnalyzer:

    def __init__(self, file_name, fitting_func, p0):
        self.x, self.y = DataAnalyzer.read_csv(file_name)
        self.fitting_func = fitting_func
        self.fitting_params = p0

    def fit(self):
        self.params, _ = curve_fit(self.fitting_func, self.x, self.y, p0=[1, 1, 0])

    def calculate_fitting_accuracy(self):
        residuals =  self.y - self.fitting_func(self.x, *self.params)
        rss = np.sum(residuals**2) 
        tss = np.sum((self.y-np.mean(self.y))**2)
        r_squared = 1 - (rss / tss)

        return r_squared
    
    def print_fitting_accuracy(self, indicator_name, value):
        print(f"-"*30)
        print(f"{indicator_name}: {value:.3}")
        print(f"-"*30)

    def evaluate_fitting_accuracy(self):
        r_squared = self.calculate_fitting_accuracy()
        self.print_fitting_accuracy("決定係数", r_squared)

    def execute(self):
        self.fit()
        self.evaluate_fitting_accuracy()

        # プロットする
        plt.scatter(self.x, self.y, label='observed data')
        plt.plot(self.x, self.fitting_func(self.x, *self.params), label=f'y={self.params[0]:.3}sin({self.params[1]:.3}x+{self.params[2]:.3})', color='red')
        plt.xlabel('x')
        plt.ylabel('y')
        plt.title('Fitting y = sin(x) to observed data')
        plt.legend()

    def read_csv(file_name):
        # CSVファイルを読み込む
        csv_file_path = file_name
        data = pd.read_csv(csv_file_path)

        # データを用意する
        x = data['x']
        y = data['y']

        return x, y

__init__という部分が状態の定義で、その後の一連のdef~が動作の定義です。
このクラスに次のようにデータを注入して、インスタンスを作成します。

def fitting_func(x, a, b, c):
    return a * np.sin(b * x + c)

data_analyzer = DataAnalyzer(
    file_name='sin_data_with_noise.csv',
    fitting_func=fitting_func,
    p0=[1, 1, 0])

data_analyzerはいろんな行動を行うことができます。
今回やりたかった一連の流れを実行することもできるし、

data_analyzer.execute()

結果
------------------------------
決定係数: 0.938
------------------------------

フィッティング精度の評価だけやることもできます。

data_analyzer.fit()
data_analyzer.evaluate_fitting_accuracy()

結果
------------------------------
決定係数: 0.938
------------------------------

インスタンス作成時に別のCSVファイルと別のフィッティング関数を用意してあげれば、一切のコピペなく、これまでの一連の流れを行うこともできます。

def fitting_func(x, a, b, c):
    return a * np.cos(b * x + c)

data_analyzer = DataAnalyzer(
    file_name='cos_data_with_noise.csv',
    fitting_func=fitting_func,
    p0=[1, 1, 0])

data_analyzer.execute()

結果
------------------------------
決定係数: 0.928
------------------------------

今回のプログラムについては、あまり恩恵を感じないかもしれないですが、クラスを使うことでプログラムの再利用性が高まることが、何となくでもわかっていただけたでしょうか。

ただし、今回作成したData_Analyzerクラスは、データの読み込みやらフィッティングやら精度評価やらプロットやら、いろんなことをやりすぎています。
例えばこの先、CSVファイルだけでなくtxtファイルやexcelファイル、その他の形式のファイルなど、さまざまな形式のファイルを読み込まなければならないかもしれません。そんな時は、データの読み込みと整形をやる別のクラスを定義するチャンスかもしれません。
また、精度評価を、もっといろんな指標でやりたくなったら、精度評価をする別のクラスを定義するときかもしれません。

このように、クラスという概念は、「プログラムの中で一定の役割を果たすもの」を定義してそれらの相互作用によってプログラムを記述していく手法に繋がっていきます。クラスをしっかり理解することが、

• 可読性
• 可変性
• 再利用性

を満たすプログラムを作成していく第一歩になると思います。

おわりに

あまりに長くなり、ここを読んでいる人がいるのか怪しいですが、今回は、段階的なプログラム作成によって、
「動くプログラムをどのように組み立てていくか」を解説し、また、関数やクラスという概念を紹介しつつ、「プログラムの再利用性」について少しでも意識を傾けよう、という話をしました。

「車輪の再発明」という言葉があるように、ソフトウェア業界では一度作られた完成度の高いものをわざわざ自分で作るなという風潮があり、一方、物理学徒には、「1から理解したい」という欲求が強い人が多いと思います。物理学徒からすれば、「中身を理解する」知的欲求と、「便利なものをブラックボックスのまま利用する」ことのいい塩梅を取りながらやるのは慣れないかもしれませんが、ともすれば良い刺激ともなるので、物理学徒の皆さん、是非プログラミングを楽しみましょう!