Использование байесовских априорных правил для штрафа за подгонку кривой ⇐ Python

[Anonymous]

Я пишу программу обработки данных для масс-спектрометрии анализа остаточных газов. В частности, я пытаюсь смоделировать потребление газа посредством ионизации как экспоненциальную функцию затухания:
y = a * exp(-p*t) + b

где:

• t — независимая переменная, время в секундах после уравновешивания газа t=0
• y — зависимая переменная, интенсивность газа в амперах.
• a, b и p — подходящие параметры.

Целью этой подгонки кривой является определение интенсивности y_0 в точке пересечения t=0. Это количество газа, которое мы измерили бы, если бы газ не нуждался в уравновешивании.
Это хорошо работает для большинства анализов, однако иногда это происходит. имеет тенденцию «переопределять» первые несколько данных и давать нереально высокие показатели потребления:


40 amu y: [1.00801174e-11 8.60445782e-12 8.74340722e-12 9.63923122e-12
8.77654502e-12 8.83196162e-12 9.44882502e-12 9.54364002e-12
8.68107792e-12 9.19894162e-12 9.26220982e-12 9.30683432e-12]
40 amu y_err: [3.55155742e-14 3.03603530e-14 3.08456363e-14 3.39750319e-14
3.09613755e-14 3.11549311e-14 3.33097888e-14 3.36410485e-14
3.06279460e-14 3.24368170e-14 3.26578373e-14 3.28137314e-14]
40 amu t: [13.489 19.117 24.829 30.433 35.939 41.645 47.253 52.883 58.585 64.292
69.698 75.408]

Обратите внимание, что прогнозируемый перехват на 11 порядков превышает любые фактические данные.
Одно из предлагаемых решений этой проблемы (которое не было выполнено) работа) заключалась в том, чтобы предоставить процедуре подгонки (curve_fit или lmfit) следующие первоначальные предположения:


Хотя это не работает, как следует из первоначальных предположений, процедура подбора кривой все равно убегает с значение p — это первоначальное предположение о p всегда гораздо более разумно, чем то, которое возвращается процедурой «оптимизации».
Я изучал байесовский метод подходы и хотели бы использовать эти первоначальные предположения и их отклонения как априорные знания для наказания совпадений, которые слишком сильно отклоняются от этих предположений - или, по крайней мере, от первоначального предположения p.
Проблема в том, что я совершенно новичок в этом и совершенно не понимаю, как это будет реализовано. Будем очень признательны за любые рекомендации.
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# exponential decay
def model(t, a, b, p):
return a * np.exp(p * t) + b

# a must be positive; p must be negative
def set_boundaries():
return [0, -np.inf, -np.inf], [np.inf, np.inf, 0]

# generate the sum term S and the first left-hand side matrix of the linear system lhs1_p
def initial_guess_p_S_and_lhs1_p(t, y):
S = [0]
for i in range(1, len(t)):
S.append(S[i-1] + 0.5 * (t - t) * (y + y))
S = np.array(S)

lhs1_p = [
[np.sum(S**2), np.sum(S*t), np.sum(S)],
[np.sum(S*t), np.sum(t**2), np.sum(t)],
[np.sum(S), np.sum(t), len(t)]
]

return S, lhs1_p

# linearize the problem and generate an initial guess of p
def initial_guess_p(t, y):
S, lhs1_p = initial_guess_p_S_and_lhs1_p(t, y)
lhs2 = [
np.sum(S*y),
np.sum(t*y),
np.sum(y)
]

init_p, _, _ = np.linalg.solve(lhs1_p, lhs2)

return init_p

# generate the left-hand side matrix of the linear system lhs1_ab
def initial_guess_ab_lhs1_ab(t, init_p):
return [
[len(t), np.sum(np.exp(init_p*t))],
[np.sum(np.exp(init_p*t)), np.sum(np.exp(2*init_p*t))]
]

# generate initial guesses of a and b
def initial_guesses_ab(t, y, init_p):
lhs1_ab = initial_guess_ab_lhs1_ab(t, init_p)

lhs2 = [
np.sum(y),
np.sum(y*np.exp(init_p*t))
]

init_b, init_a = np.linalg.solve(lhs1_ab, lhs2)

return init_a, init_b

# generate the residuals of the model
def initial_residuals(t, y, a, b, p):
return y - model(t, a, b, p)

# generate the initial guess of the variance of the residuals
def init_p_error(t, y, init_p, initial_resids):
S, lhs1_p = initial_guess_p_S_and_lhs1_p(t, y)
init_p_variance, _, _ = np.linalg.inv(lhs1_p) * np.sum(initial_resids**2)

return init_p_variance

# generate the initial guess of the variance of the residuals
def init_ab_error(t, init_p, initial_resids):
lhs1_ab = initial_guess_ab_lhs1_ab(t, init_p)
init_b_variance, init_a_variance = np.linalg.inv(lhs1_ab) * np.sum(initial_resids**2)

return init_a_variance, init_b_variance

# Add the main function to test the above functions
def main():
t = np.linspace(0, 10, 100)
a, b, p = 2.5, 1.0, -0.5
y = model(t, a, b, p) + 0.1 * np.random.normal(size=t.size)

init_p = initial_guess_p(t, y)
init_a, init_b = initial_guesses_ab(t, y, init_p)

print(f"Initial guesses: a={init_a}, b={init_b}, p={init_p}")

# Fit the model with initial guesses
popt, _ = curve_fit(model, t, y, p0=[init_a, init_b, init_p])

# Plotting
plt.figure()
plt.scatter(t, y, label='Data')
plt.plot(t, model(t, *popt), label='Fit', color='red')
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.show()

if __name__ == "__main__":
main()


Подробнее здесь: https://stackoverflow.com/questions/792 ... ve-fitting