検量線を引き、吸光度から濃度を求める¶
この作業がかなり多く、パターン化しているので、まとめておく。なお、可読性のために変数名は長めにしている。
How to use:¶
$ python Calibration_curve.py --Concentration 0,0,0.10,0.25,0.50,1.00,1.50,2.00,2.50 \
--Absorbance 0.048,0.042,0.045,0.045,0.076,0.103,0.151,0.202,0.262 \
--unit mM \
--wavelength 540
import numpy as np
import matplotlib.pyplot as plt
実験データ¶
| 濃度(mg/mL) | 0 | 0.01 | 0.02 | 0.05 | 0.10 | 0.20 | 0.40 | 0.60 | 0.80 | 1.00 | 2.00 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 吸光度($A_{280}$) | 0 | 0.004 | 0.010 | 0.000 | 0.043 | 0.067 | 0.159 | 0.288 | 0.382 | 0.473 | 0.962 |
前処理¶
unit = "mg/mL"
wavelength = "$A_{280}$"
# 上の表をそのままベタ貼り。
c_str = "濃度(mg/mL) 0 0.01 0.02 0.05 0.10 0.20 0.40 0.60 0.80 1.00 2.00"
a_str = "吸光度( A280 ) 0 0.004 0.010 0.000 0.043 0.067 0.159 0.288 0.382 0.473 0.962"
# データ部分だけを取得する。
c_ = [float(data) for data in c_str.split("\t")[1:]]
a_ = [float(data) for data in a_str.split("\t")[1:]]
# ここまでのデータ処理（表がないことがほとんどなので、これを最初から作っても良い）
print("Concentration: {}".format(c_))
print("Absorbance : {}".format(a_))
# sklearn が扱えるデータの型に直す。
c = np.array(c_).reshape(-1, 1)
a = np.array(a_).reshape(-1, 1)
学習¶
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error as mse
lr = LinearRegression()
# 学習
lr.fit(c,a)
coef = lr.coef_[0][0] # 係数
intercept = lr.intercept_[0] # y切片
# 検量線
Calibration_curve = "$y={:.3f}x{}{:.3f}$".format(coef, "+" if intercept > 0 else "", intercept)
# 決定係数
Coefficient_of_determination = "$R^2$ = {:.3f}".format(lr.score(c,a))
# 平均二乗誤差
Root_Mean_Square_Error = "RMSE = {:.3f}".format(mse(a, lr.predict(c)) ** (1/2))
# プロット
plt.plot(c, lr.predict(c), label="Calibration curve: "+Calibration_curve)
plt.scatter(c,a, label="{}\n{}".format(Coefficient_of_determination,
Root_Mean_Square_Error))
plt.xlabel("Concentration [{}]".format(unit))
plt.ylabel("Absorbance ({})".format(wavelength))
plt.title("The relationship between Concentration and Absorbance")
plt.legend()
plt.show()
外れ値を除く¶
上記の結果では、濃度 0.05 mg/mL時に吸光度が 0 となっているのは、明らかにヒューマンエラーだと分かる。そこで、この値を外れ値だとみなし、残りのデータでフィッティングを行う。
# インデックスを利用して
Outlier_Concentration = 0.05
[c for c,a in zip(c_,a_) if c != Outlier_Concentration]
c_data = np.array([c for c,a in zip(c_,a_) if c != Outlier_Concentration]).reshape(-1, 1)
c_outlier = np.array([c for c,a in zip(c_,a_) if c == Outlier_Concentration]).reshape(-1, 1)
a_data = np.array([a for c,a in zip(c_,a_) if c != Outlier_Concentration]).reshape(-1, 1)
a_outlier = np.array([a for c,a in zip(c_,a_) if c == Outlier_Concentration]).reshape(-1, 1)
lr_2 = LinearRegression()
lr_2.fit(c_data, a_data)
coef_2 = lr_2.coef_[0][0] # 係数
intercept_2 = lr_2.intercept_[0] # y切片
Calibration_curve_2 = "$y={:.3f}x{}{:.3f}$".format(coef_2, "+" if intercept_2 > 0 else "", intercept_2)
Coefficient_of_determination_2 = "$R^2$ = {:.3f}".format(lr_2.score(c_data,a_data))
Root_Mean_Square_Error_2 = "RMSE = {:.3f}".format(mse(a_data, lr_2.predict(c_data)) ** (1/2))
# プロット
plt.plot(c_data, lr_2.predict(c_data), label="Calibration curve: "+Calibration_curve_2)
plt.scatter(c_data,a_data, label="{}\n{}".format(Coefficient_of_determination_2,
Root_Mean_Square_Error_2))
plt.scatter(c_outlier, a_outlier, label="Outlier")
plt.xlabel("Concentration [{}]".format(unit))
plt.ylabel("Absorbance ({})".format(wavelength))
plt.title("The relationship between Concentration and Absorbance")
plt.legend()
plt.show()
検量線を関数化する¶
上のグラフとは異なり、吸光度から濃度を求める関数を求めたいので、上の検量線の逆関数を求める。
A2C = lambda y: (y-intercept)/coef
あとは、この関数を未知試料の吸光度に適応すれば未知試料の濃度が求まる。
※ここからは、実際のレポートに書いた考察をつらつらと書いています。間違いがあれば全体を通してですが、指摘していただけると嬉しいです。
重み付き最小二乗誤差¶
上記で記した回帰直線は、全てのデータ点 $y_i$ が $\hat{y_i}$ を平均とした正規分布に従って出現するという仮定の下で誤差を最小にした直線であり、それがゆえに決定係数 $R^2$ が意味を成す。しかし、実際のデータでは、どの仮定は一般に成り立たない。その理由の簡単な説明は、「使っているピペットマンが違うから」である。
どのように希釈系列を作ったかにもよるが、本実験では希釈系列の総量は揃えている。従って、希釈試料を調整する際に誤差に大きな影響を及ぼすのは水の量ではなく試料の量である。
ここで、ピペットマンにはいくつかの種類があり、それぞれ測り取れる最大容量が異なる。測定時の誤差はピペットマンの最大容量に比例するはずなので、以上の理論により「希釈系列の濃度が高いほど誤差が大きい」と結論づけられる。
厳密には、それぞれの希釈系列を作製する際にどのピペットマンを使ったか、段階希釈をしたか、などの情報が必要になるが、ここでは簡単のために「希釈濃度が誤差の大きさに比例する」と仮定する。
定式化¶
以上を踏まえると、希釈系列の濃度 $\{x_i\}$ 吸光度 $\{y_i\}$ のデータが得られた時、ナイーブな最小二乗では $$E=\frac{1}{n}\sum_{i}\left\{y_{i}-\left(a x_{i}+b\right)\right\}^{2}$$ を最小化していたが、今回の最小二乗では $$E^{\prime}=\frac{1}{n}\sum_{i}\left\{\frac{y_{i}-\left(a x_{i}+b\right)}{x_{i}}\right\}^{2}$$ を最小化する必要がある。
通常の誤差関数と同様に $a,b$ で偏微分をして $0$ とおけば解析的に解くことができ、結果は以下のようになる。
$$\left\{\begin{array}{l}{a=\frac{n}{n^{2}-1}\left\{\sum_{i} \frac{y_{i}}{x_{i}}-\frac{1}{n} \sum_{i} \frac{1}{x_{i}} \sum_{i} y_{i}\right\}} \\ {b=\frac{n}{n^{2}-1}\left\{\sum_{i} y_{i}-\frac{1}{n} \sum_{i} \frac{y_{i}}{x_{i}} \sum_{i} x_{i}\right\}}\end{array}\right.$$なお、この際決定係数を使って回帰直線の評価を行うことができないが、誤差関数の値をそのまま比較すれば良いことになる。
ブランクの取り扱い¶
※ この際正しく数値計算するためには、$x_i=0$ となる値（ブランク）の処理が非常に大事になる。
ここでは、「ブランクには誤差が存在しない」という大胆な仮定をおき、回帰直線がブランクを必ず通るようにすることで、この問題に対処する。
ブランクのデータを $(x_{\mathrm{blank}}, y_{\mathrm{blank}})$ とすると、この点を通る回帰直線は以下のようにかける。
$$y-y_{\mathrm{blank}} = m\left(x-x_{\mathrm{blank}}\right)$$すると、誤差関数は以下のようにかける。 $$E_{\mathrm{blank}}=\frac{1}{n}\sum_{i}\left\{\frac{y_{i}-\left\{m\left(x_i-x_{\mathrm{blank}}\right) + y_{\mathrm{blank}}\right\}}{x_i}\right\}^{2}$$
この誤差を最小化する $m$ は解析的に解くことができ、以下で求まる。 $$m = \frac{\sum_i\frac{\left(y_i-y_{\mathrm{blank}}\right)\left(x_i-x_{\mathrm{blank}}\right)}{x_i^2}} {\sum_i\frac{\left(x_i-x_{\mathrm{blank}}\right)^2}{x_i^2}}$$
実装¶
「実装」と言う程ではないが、先ほどのデータに適用してみる。
x = c[1:]
y = a[1:]
x_blank = c[0]
y_blank = a[0]
m = sum((y-y_blank)*(x-x_blank)/x**2) / sum((x-x_blank)**2/x**2)
coef_3 = m[0] # 係数
intercept_3 = (coef_3 * (0-x_blank) + y_blank)[0] # y切片
Calibration_curve_3 = "$y={:.3f}x{}{:.3f}$".format(coef_3, "+" if intercept_3 >= 0 else "", intercept_3)
func3 = lambda x:coef_3*(x-x_blank) + y_blank
plt.plot(c, lr.predict(c), label="Calibration curve: "+Calibration_curve)
plt.plot(x, np.vectorize(func3)(x), label="Calibration curve (weighted): "+ Calibration_curve_3)
plt.scatter(c,a,label="data")
plt.xlabel("Concentration [{}]".format(unit))
plt.ylabel("Absorbance ({})".format(wavelength))
plt.title("The relationship between Concentration and Absorbance")
plt.legend()
plt.show()
データ点の偏りを均等に捉えている気がする。（が、この例だと元の回帰直線の方が綺麗に近似できているように見える。）