התאמת מודלים והתפלגויות בעזרת SciPy

במחברת זו נראה איך SciPy עוזרת לנו לעבור מנקודות מדידה למודל כמותי. נתחיל בעקומת כיול פשוטה, נבדוק את השגיאות בעזרת היסטוגרמה, ואז נבצע התאמה לא־לינארית של דעיכה קינטית.

המטרה אינה ללמוד סטטיסטיקה מתקדמת, אלא לראות תבנית עבודה שימושית: מגדירים מודל, משווים אותו לנתונים, ומסתכלים על מה שנשאר.

🤖 Ask PyProf.

🗑 −

▶ שלח

⏳ ממתין...

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

דוגמה 1: עקומת כיול

נניח שמדדנו בליעה של תמיסות בריכוזים ידועים. לפי חוק בר־למברט, בתחום מתאים, הבליעה אמורה להיות בקירוב לינארית בריכוז:

\[A = a c + b\]

כאן \(c\) הוא הריכוז, \(A\) היא הבליעה, \(a\)הוא השיפוע, ו־ \(b\) הוא חיתוך עם הציר. בפועל נקבל גם רעש מדידה.

ראשית, נחולל נתונים סינתתיים האמורים לדמות תוצאה של מדידות ניסיוניות.

rng = np.random.default_rng(42)

concentration = np.linspace(0.0, 2.0, 2000) # mmol/L
true_slope = 0.75
true_intercept = 0.04
noise = rng.normal(loc=0.0, scale=0.07, size=concentration.size)

absorbance = true_slope * concentration + true_intercept + noise

שנית, נשרטט את הנתונים שלנו.

fig, ax = plt.subplots()
ax.scatter(concentration, absorbance, s=0.3)
ax.set_xlabel('concentration (mmol/L)')
ax.set_ylabel('absorbance')
ax.set_title('Calibration data')
plt.show()

כדי לבצע התאמה בעזרת curve_fit, קודם מגדירים פונקציה של המודל. הארגומנט הראשון הוא המשתנה הבלתי־תלוי, ואחריו הפרמטרים שאותם רוצים למצוא.

def linear_model(c, slope, intercept):
    return slope * c + intercept

params, covariance = curve_fit(linear_model, concentration, absorbance)
slope, intercept = params

print(f'slope = {slope:.3f}')
print(f'intercept = {intercept:.3f}')

slope = 0.745
intercept = 0.041

c_fit = np.linspace(0, 2.0, 200)
A_fit = linear_model(c_fit, slope, intercept)

fig, ax = plt.subplots()
ax.scatter(concentration, absorbance, label='measurements',s=0.3)
ax.plot(c_fit, A_fit, label='fit', color='black')
ax.set_xlabel('concentration (mmol/L)')
ax.set_ylabel('absorbance')
ax.set_title('Linear calibration fit')
ax.legend()
plt.show()

כעת אפשר להשתמש בעקומת הכיול כדי להעריך ריכוז של דוגמה לא־ידועה. אם הבליעה של הדוגמה היא \(A_\mathrm{unknown}\), אז לפי המודל הלינארי:

\[c_\mathrm{unknown} = \frac{A_\mathrm{unknown} - b}{a}\]

A_unknown = 0.62
c_unknown = (A_unknown - intercept) / slope

print(f'unknown concentration = {c_unknown:.3f} mmol/L')

unknown concentration = 0.777 mmol/L

שאריות והיסטוגרמה

השארית היא ההפרש בין המדידה לבין ערך המודל באותה נקודה. אם המודל סביר והרעש אקראי, השאריות לא אמורות להראות מבנה ברור. כמובן, כיוון שבמקרה שלנו הוספנו רעש גאוסי מלאכותי, זה מה שאנחנו מצפים לראות.

predicted = linear_model(concentration, slope, intercept)
residuals = absorbance - predicted

fig, ax = plt.subplots()
ax.scatter(concentration, residuals, s=0.3)
ax.axhline(0, linestyle='--', color='black')
ax.set_xlabel('concentration (mmol/L)')
ax.set_ylabel('residual')
ax.set_title('Residuals')
plt.show()

במקום להביט בעין, אפשר לעשות בדיקה שיטטית יותר באמצעות היסטוגרמה. יש פונקציה ב־matplotlib שעושה את זה עבורנו.

mu, sigma = norm.fit(residuals)
x = np.linspace(residuals.min() - 0.02, residuals.max() + 0.02, 200)

fig, ax = plt.subplots()
ax.hist(residuals, bins=50, density=True, alpha=0.6)
ax.plot(x, norm.pdf(x, mu, sigma), label='normal fit')
ax.set_xlabel('residual')
ax.set_ylabel('density')
ax.set_title('Histogram of residuals')
ax.legend()
plt.show()

print(f'mean residual = {mu:.4f}')
print(f'std residual = {sigma:.4f}')

mean residual = -0.0000
std residual = 0.0701

דוגמה 2: התאמה לא־לינארית לדעיכה קינטית

כעת נשתמש באותה פקודה, אבל עם מודל לא־לינארי. נניח שמדדנו בליעה של חומר שמגיב לפי קינטיקה מסדר ראשון. אם הבליעה פרופורציונית לריכוז, נקבל:

\[A(t) = A_\infty + (A_0 - A_\infty)e^{-kt}\]

כאן \(k\) הוא קבוע הקצב, ו־\(A_\infty\) הוא ערך הבליעה לאחר זמן ארוך.

נשתמש כאן שוב ב־np.exp, הפונקציה של NumPy שמחשבת \(e^x\) על מספרים או על מערכים. במקרה הזה np.exp(-k * t) נותנת את גורם הדעיכה לכל הזמנים במערך t.

time = np.linspace(0, 10, 1000) # min
true_A0 = 1.20
true_Ainf = 0.18
true_k = 0.32

def decay_model(t, A0, Ainf, k):
    return Ainf + (A0 - Ainf) * np.exp(-k * t)

abs_decay = decay_model(time, true_A0, true_Ainf, true_k)
abs_decay += rng.normal(0.0, 0.025, size=time.size)

guess = [1.0, 0.1, 0.2]
params, covariance = curve_fit(decay_model, time, abs_decay, p0=guess)
A0_fit, Ainf_fit, k_fit = params

print(f'A0 = {A0_fit:.3f}')
print(f'Ainf = {Ainf_fit:.3f}')
print(f'k = {k_fit:.3f} 1/min')

A0 = 1.200
Ainf = 0.182
k = 0.320 1/min

t_fit = np.linspace(0, 10, 300)
A_fit = decay_model(t_fit, A0_fit, Ainf_fit, k_fit)

fig, ax = plt.subplots()
ax.scatter(time, abs_decay, label='measurements', s=0.3)
ax.plot(t_fit, A_fit, label='nonlinear fit', color='black')
ax.set_xlabel('time (min)')
ax.set_ylabel('absorbance')
ax.set_title('First-order kinetic decay')
ax.legend()
plt.show()

שימו לב שהתהליך המחשבתי דומה מאוד להתאמה הלינארית: מגדירים מודל, נותנים ל־SciPy את הנתונים, ומקבלים פרמטרים. ההבדל הוא שהקשר בין הפרמטרים לבין הנתונים כבר אינו קו ישר.

תרגילים קצרים

1. הגדילו את רעש המדידה בעקומת הכיול ובדקו מה קורה לשיפוע המותאם.

2. שנו את \(k\) בדעיכה הקינטית וראו איך משתנה הגרף.

3. נסו לבצע התאמה לדעיכה בלי לתת ניחוש התחלתי p0. האם התוצאה משתנה?