使用贝叶斯先验来惩罚曲线拟合

问题描述 投票：0回答：1

我正在编写一个用于残余气体分析质谱的数据缩减程序。具体来说，我正在尝试通过电离将气体消耗建模为指数衰减函数：

y = a * exp(-p*t) + b

地点：

  • t
    是自变量，气体平衡后的时间（以秒为单位） t=0
  • y
    是因变量，气体强度（以安培为单位）
  • a
    、
    b
    和 
    p
    是拟合参数

此曲线拟合的目标是确定截距 t=0 处的强度 

y_0
。这代表了如果气体不需要平衡的话我们会测量的气体量。 这对于大多数分析都很有效，但是，有时它有“过度拟合”前几个数据并产生不切实际的高消耗率的倾向：

enter image description here 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）：

enter image description here 虽然这不能作为初始猜测——曲线拟合例程仍然会偏离

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()



编辑：添加更多示例


He 4390 21-TC-54-4 H11 RE
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]



He 4399: CB005
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
 及其错误作为贝叶斯先验。
    

	




python

curve-fitting

bayesian


1个回答




    
    
    

        

            

                

                    
0
投票

                    
                

            

        

        
curve_fit
 在给定蹩脚数据的情况下尽力而为的结果。
因此，解决方案是“不是”尝试用贝叶斯方法“驯化”指数衰减模型，而是完全避免拟合这种嘈杂的趋势。这是通过使用 Aikake 信息标准计算两个竞争模型来实现的：指数衰减模型和总体平均值。程序使用 AIC 最低的方法。


在所有示例情况下，将指数 AIC 与平均 AIC 进行比较可以得出“平均值”在所有情况下都是更好的方法。

    
	



    

    


  

  

    
最新问题


  





  

    © www.soinside.com 2019 - 2024. All rights reserved.