Использование байесовских априорных правил для штрафа за подгонку кривой ⇐ 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[i] - t[i - 1]) * (y[i] + y[i - 1]))
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()

РЕДАКТИРОВАНИЕ: добавление дополнительных примеров
[img]https:/ /i.sstatic.net/5KTiJFHO.png[/img]

Код: Выделить всё

40 amu p init: -0.14367518985061087 p: -0.6565324927370925
40 amu y: [1.03659297e-11 9.01981636e-12 9.21698536e-12 8.42299856e-12  8.99095476e-12 9.32831176e-12 8.66239446e-12 8.94233696e-12  9.66985026e-12 8.93559296e-12 9.08987596e-12 8.93206396e-12]
40 amu y_err: [4.58060424e-14 3.98911199e-14 4.07573504e-14 3.72694684e-14 3.97643259e-14 4.12464696e-14 3.83209953e-14 3.95507423e-14 4.27471417e-14 3.95211154e-14 4.01989091e-14 3.95056122e-14]
40 amu t: [ 8.908 14.511 20.117 25.833 31.434 37.146 42.648 48.249 53.954 59.586
65.211 70.822]

Код: Выделить всё

40 amu p init: -0.15494649546677705 p: -1.7804944874902089
40 amu y: [1.02342173e-11 9.13542299e-12 8.71434679e-12 9.30896839e-12 9.67921739e-12 8.81455689e-12 9.01517339e-12 9.32982869e-12  9.07950499e-12 9.10221369e-12 9.13479289e-12 9.74699459e-12]
40 amu y_err: [4.52082398e-14 4.03777891e-14 3.85269625e-14 4.11406522e-14 4.27682641e-14 3.89674163e-14 3.98492179e-14 4.12323508e-14 4.01319933e-14 4.02318124e-14 4.03750194e-14 4.30662243e-14]
40 amu t: [ 9.808 15.54  21.056 26.757 32.365 37.967 43.603 49.221 54.934 60.453
66.158 71.669]

Обратите внимание, что в обоих случаях начальное предположение p меньше и, следовательно, более разумно, чем любое значение, полученное Curve_fit. Поэтому я хочу использовать это первоначальное предположение p и его ошибку в качестве байесовского априора.

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