(学习笔记）python 分区间统计，加权最小二乘法拟合数据：xlrd，sicpy.optimize.leastsq

最新推荐文章于 2022-09-07 13:29:49 发布

greatofdream

最新推荐文章于 2022-09-07 13:29:49 发布

阅读量3k

点赞数 2

分类专栏： python 文章标签： python

本文链接：https://blog.csdn.net/zaq15csdn/article/details/85163808

版权

python 专栏收录该内容

12 篇文章 2 订阅

订阅专栏

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

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

$\chi^2=\sum{(n_i-m_i)^2/\sigma}$

计算，其中分母项都是 $\sigma^2$ ,所以可以忽略，直接极小化。这个公式一般适用于很多情况，因为噪声大部分情况是和采样点无关的。

对于通常的计数观测，就是每个bin里统计有多少个事例，每个bin里服从泊松分布，所以对应的方差就是观测的值 $\sigma^2=n_i$ ，因此最小二乘法下面的分母可以替换成对应的方差值。这里的误差实际上是统计误差，而且没有考虑其它的本底。

最小二乘法 $\chi^2=\sum{(n_i-m_i)^2/n_i}$

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)
#绘制拟合曲线