3

【最小二乘估计】scipy.optimize.leastsq

 3 years ago
source link: https://www.guofei.site/2017/06/06/scipyleastsq.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
neoserver,ios ssh client

【最小二乘估计】scipy.optimize.leastsq

2017年06月06日

Author: Guofei

文章归类: 5-6-最优化 ，文章编号: 7301

版权声明：本文作者是郭飞。转载随意，但需要标明原文链接，并通知本人
原文链接：https://www.guofei.site/2017/06/06/scipyleastsq.html

Edit

常见的曲线拟合

趋势模型里写了趋势模型中常用的10种曲线

这里是曲线拟合种常用的曲线

  1. Polynomial Models
    y=∑i=1npixiy=∑i=1npixi
  2. Exponential Models
    y=aebx,y=aebx+cedxy=aebx,y=aebx+cedx
  3. Fourier Series y=a+∑i=1naicos(nwx)+bisin(nwx)y=a+∑i=1naicos⁡(nwx)+bisin⁡(nwx)
  4. Gaussian Models
    y=∑i=1naiexp(−(x−bici)2)y=∑i=1naiexp⁡(−(x−bici)2)
  5. Power Series
    y=axb,y=a+bxcy=axb,y=a+bxc
  6. Rational Polynomials
    y=∑i=1npixi∑i=1nqixiy=∑i=1npixi∑i=1nqixi
  7. Sum of Sines Models y=a+∑i=1naisin(nwx)y=a+∑i=1naisin⁡(nwx)
  8. Weibull Distributions y=abxb−1e−axby=abxb−1e−axb
    widely used in reliability and life (failure rate) data analysis

leastsq

最小二乘估计原理是这样的：

y=f(x,θ)+εy=f(x,θ)+ε,
其中εε独立同分布。

θ=argmin∑(yi−f(xi,θ))2θ=arg⁡min∑(yi−f(xi,θ))2

非线性最小二乘法中，SST=SSR+SSE不再成立，但仍然可以定义R_squared=1-SSE/SST

leastsq可以用来做最小二乘估计，可以在线性拟合和非线性拟合中使用。

线性拟合案例

step1:生成模拟的源数据

import numpy as np
from scipy.optimize import leastsq

# 生成模拟数据
X=np.linspace(10,40,1000)
Y=0.8*X+2.1+np.random.normal(loc=0,scale=1,size=[1,1000])[0]

step2：定义残差并寻优

def residuals(p):
    k,b=p
    return Y-(k*X+b)
r=leastsq(residuals,[1,0])#[1,0]为初始值
k,b=r[0]
k,b

step3:计算误差

def S(k,b):
    error=np.zeros(k.shape)
    for x,y in zip(X,Y):
        error+=(y-(k*x+b))**2
    return error
S(k,b)

综合案例

step1：定义几个目标函数，准备一一检验其拟合效果

import numpy as np

def test_func1(x,p):
    a,b=p
    return a*x+b

def test_func2(x,p):
    A,k,theta=p#这是一个编程技巧，虽然冗余了一行，但是可读性极大提高
    return A*np.sin(2*np.pi*k*x+theta)

def test_func3(x,p):
    a,b,c=p
    return a*x**2+b*x+c

step2:定义目标函数和残差函数

p_true=[0.4,-2,0.9]#真实值
# 目标函数
def obj_func(x,p):
    return test_func3(x,p)

# 残差
def residuals(p,y,x):
    return y-obj_func(x,p)

step3:生成模拟数据

X=np.linspace(0,10,1000)
y=obj_func(X,p_true)+np.random.randn(len(X))

step4:拟合

from scipy.optimize import leastsq
p_prior=np.ones_like(p_true)# 先验的估计，真实数据分析流程中，先预估一个接近的值。这里为了测试效果，先验设定为1

plsq=leastsq(residuals,p_prior,args=(y,X))

print(p_true)
print(plsq)

step5:画图

import matplotlib.pyplot as plt
plt.plot(X,y)
plt.plot(X,obj_func(X,plsq[0]))
plt.show()

下面这3个目标函数拟合效果的展示：
fun1.png
fun2.png
fun3.png

您的支持将鼓励我继续创作！

report

Recommend

About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK