多元加权最小二乘估计python代码 python 加权最小二乘

自己用python3.x处理数据遇到的问题，在这里记录分享一下。

最小二乘法研究的问题是y=Ax+n，其中y是观测值，x是采样点，n是噪声，A是需要拟合的系数矩阵，通常我们认为噪声是白噪声，所以n服从正态分布N~(0,

)，那么我们在计算最小二乘法时对

计算，其中分母项都是

,所以可以忽略，直接极小化。这个公式一般适用于很多情况，因为噪声大部分情况是和采样点无关的。对于通常的计数观测，就是每个bin里统计有多少个事例，每个bin里服从泊松分布，所以对应的方差就是观测的值

，因此最小二乘法下面的分母可以替换成对应的方差值。这里的误差实际上是统计误差，而且没有考虑其它的本底。最小二乘法

1.首先从excel读取数据 ，使用python的包xlrd

import xlrd
#导入读取excel包

f=xlrd.open_workbook('jisuan.xlsx')
#创建一个文件指针，此处文件名为jisuan.xlsx

sheet=f.sheet_by_index(1)
#读取excel的第2张表

tdata=sheet.col_values(1)
#读取该表的第一列

2.使用python的包pandas进行分组统计

此处只找到了一个cut函数，可能会有更好的方法。

matlab中有一个histcounts方法还挺好用的，此处就不详细叙述了

tmax=max(tdata)
tmin=min(tdata)
#找到最大最小值0.3，8.8

tbin=[0.3+i*0.4 for i in range(23)]
#这个区间是自己算了一下，用0.4作为区间间隔，0.3开始

fpinshu=pd.cut(tdata,tbin,right=False)
#用pandas进行区间统计

fnum=fpinshu.value_counts()
#取出频数，很可惜的是fnum不是list，但是可以索引

3.使用python的包scipy进行拟合

import numpy as np
from scipy.optimize import leastsq

def func(p,x):
    A,tau,gamma=p
    return A*exp(-x/tau)+gamma
#我要拟合的函数模型为y=A*exp(-t/tau)+gamma

def error(p,x,y):
    temp=func(p,x)
    return (temp-y)/np.sqrt(temp)
#定义需要的误差函数，需要在leastsq中用到.注意这里用的sqrt和exp都是numpy的函数，可以支持数组运算如果报错可以考虑是不是import了math中的sqrt和exp。

start=[1,2.2,1]
#定义三个参数的初值

para=leastsq(error,start,args=(tx[0:15],y[0:15]))
#tx是区间横坐标的中点，y是计数，因为拟合时出现了被除数过小的问题，直接把后面的数据舍去，用前面16个点拟合。

可以参考
https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.optimize.leastsq.html

4.使用matplotlib 画图

import matplotlib.pyplot as plt

n,bins,pathces=plt.hist(t,tbin)
#此处可以直接得到频数n，也就是说上面使用pandas其实是多次一举

y=func(para[0],tx)
#得到拟合结果y

plt.plot(x,y)
#绘制拟合曲线

最后的拟合结果还是有一些问题，出现了小于0的情况，应该是数据量不够的原因。