最小二乘法

Question

我有一组数据，它遵循双高斯分布，它的方程依赖于4个参数。我想做的是：

• 从4个参数的值开始(这些参数与我从数据集中获得的双高斯参数更接近)
• 逐步增加
• 用这些参数绘制双高斯图
• 检查，使用最小二乘法，这是更好的值，绘制一个双高斯，这是最适合的双高斯，我得到的数据集。

我试过这样做，下面是我的代码：

import matplotlib.pyplot as plt
import numpy as np

def DGauss(x,I1,I2,sigma1,sigma2):
    return (I1/np.sqrt(2*np.pi*sigma1))*np.exp(-x*x/(2*sigma1*sigma1)) + (I2/np.sqrt(2*np.pi*sigma2))*np.exp(-x*x/(2*sigma2*sigma2))

def Inte(x0,xf,npoints,I1,I2,sigma1,sigma2):
     delta = (xf-x0)/npoints
     h = 0
     for i in range(npoints):
         x = x0 + i*delta
         if (i % 2) == 0:
             fact = 4.0
         else:
             fact = 2.0 

         h += fact*DGauss(x,I1,I2,sigma1,sigma2)

     return (abs(delta)/3.0)*(h + DGauss(x0,I1,I2,sigma1,sigma2) + DGauss(xf,I1,I2,sigma1,sigma2))

def IntBLM(start,end):
    delta = 1
    VarIntBLM = (1.0/3.0)*(BLM[start]+4.0*BLM[start+delta]+BLM[end])
    return VarIntBLM

Position = [9.88, 8.68, 7.48, 6.28, 5.08, 3.88, 3.28, 3.13, 3.08, 3.03, 2.98, 2.93, 2.88, 2.83, 2.78, 2.73, 2.68, 2.63, 2.58, 2.53, 2.48, 2.43, 2.38, 2.33, 2.28, 2.23, 2.18, 2.13, 2.08, 2.03, 1.98, 1.93, 1.88, 1.83, 1.78, 1.73, 1.68, 1.63, 1.58, 1.53, 1.48, 1.43, 1.38, 1.33, 1.28, 1.23, 1.18, 1.13, 1.08, 1.03, 0.98, 0.93, 0.88, 0.83, 0.78, 0.73, 0.68, 0.63, 0.58, 0.53, 0.48, 0.43, 0.38, 0.33, 0.28, 0.23, 0.18, 0.13, 0.08, 0.03]

yVal = [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.05493778e-04, 3.68936803e-04, 1.73977308e-03, 9.86279398e-03, 1.52954321e-02, 2.42624032e-02, 2.87455974e-02, 3.23848413e-02, 3.28592719e-02, 3.94523416e-02, 4.61509051e-02, 5.70161813e-02, 6.37672003e-02, 7.19426767e-02, 7.76393407e-02, 8.56568678e-02, 9.61526244e-02, 1.04328101e-01, 1.13506059e-01, 1.19940597e-01, 1.26006198e-01, 1.40933276e-01, 1.50796653e-01, 1.66514643e-01, 1.80650226e-01, 1.93889404e-01, 2.04754098e-01, 2.17940237e-01, 2.28067057e-01, 2.37930434e-01, 2.51644042e-01, 2.63511800e-01, 2.80759742e-01, 2.95686820e-01, 3.08715010e-01, 3.31184602e-01, 3.46480617e-01, 3.69846614e-01, 3.85406655e-01, 4.06188347e-01, 4.28394495e-01, 4.50020137e-01, 4.83039107e-01, 5.07460625e-01, 5.31670573e-01, 5.54879204e-01, 5.78351278e-01, 6.02561809e-01, 6.25664363e-01, 6.57048471e-01, 6.82893863e-01, 7.13327944e-01, 7.32580267e-01, 7.69608000e-01, 7.87699892e-01, 8.14072753e-01, 8.33588519e-01, 8.52102386e-01, 8.71090683e-01, 8.94562175e-01, 9.16187816e-01, 9.37602470e-01, 9.56802338e-01, 9.69197565e-01, 9.78321903e-01, 9.84861934e-01, 9.93142904e-01,]

BLM = [0., 0., 0., 0.000181, 0.000452, 0.002352, 0.000724, 0.000633, 0.001267, 0.000633, 0.006244, 0.000814, 0.002805, 0.000633, 0.002443, 0.000814,0.001629, 0.001086, 0.001448, 0.001357, 0.002443, 0.001991, 0.002805, 0.010407,0.001991, 0.000814, 0.002352, 0.002714, 0.002081, 0.002714, 0.002352, 0.00181, 0.004253, 0.001719, 0.003077, 0.001448, 0.002805, 0.002896, 0.002986, 0.0019, 0.002262, 0.002262, 0.004072, 0.004434, 0.006425, 0.002986, 0.005701, 0.002352,0.004163, 0.004615, 0.001991, 0.001538, 0.005158, 0.0019,   0.003981, 0.004615,0.004615, 0.003167, 0.005339, 0.003619, 0.005339, 0.00181,  0.002624, 0.001267, 0.0038, 0.001629, 0.002986, 0.000633, 0.001719, 0.000543]

npoints = 13 
npairs = 13


hMin = 1
h=0

for k in range(13):
    sigma1 = 0.01*(k+1) + 1
    for l in range(13):
        sigma2 = 0.01*(l+1) + 1
        for m in range(13):
            I1 = 0.1*(m+1)
            for j in range(13):
                I2 = 0.1*(j+1)
                for i in range(13):
                    startPos = 50+i 
                    endPos = 52+i 
                    startX = Position[startPos]
                    endX = Position[endPos]

                    h += (IntBLM(startPos,endPos)-Inte(startX,endX,npoints,I1,I2,sigma1,sigma2))**2

                if h < hMin:
                    hMin = h
                    I1min = I1
                    I2min = I2
                    sigma1min = sigma1
                    sigma2min = sigma2

plt.figure(1)
plt.plot(Position,yVal,'o')
plt.plot(Position,DGauss(Posi,I1min,I2min,sigma1min,sigma2min), label='Least Square')
plt.show()

我做了不同的测试，打印高斯参数的值，最小二乘法返回的值，但问题是，I1min, I2min, sigma1min, sigma2min的值总是来自循环的第一个循环，I1min的值总是等于I2min，链接sigma1min的值等于sigma2min，我知道这是不对的。有人能帮帮我吗？谢谢!

Answer 1

好的，逻辑是这样的：

1. hMin = 1，h=0
   1. h += (IntBLM(startPos,endPos)-Int(startX,endX,npoints,I1,I2,sigma1,sigma2))**2 a. 如果 h大于或等于1，则条件永不执行。 如果 h现在介于0到1之间，则条件执行1. hMin = h。 现在，h将永远不会小于hMin，因为它总是增加一个正数(一个平方值)。因此，条件不再执行，并且I1min, I2min, sigma1min, sigma2min维护它们的实例化值。

​

​