# -*- coding: utf-8 -*8
"""
Created on Wed Dec 11 21:22:24 2019
@author: Admin
"""
##最小二乘法
import numpy as np ##科学计算库
import scipy as sp ##在numpy基础上实现的部分算法库
import matplotlib.pyplot as plt ##绘图库
import math
from scipy.optimize import leastsq ##引入最小二乘法算法
##样本数据(Xi,Yi),需要转换成数组(列表)形式
Xi=np.array([10,20,30,40,50,60,70,80,90,100])
Yi=np.array([0.00233982,0.0306378,0.137365,0.363723,0.661229,0.890883,0.982852,
0.99898,0.999983,0.9999999])
Yj=np.log(np.log(1/(1-Yi)))
Xj=np.log(Xi)
#k=m
#b=-m*np.log(YITA)
##需要拟合的函数func :指定函数的形状
def func(p,x):
k,b=p
#m=k #m为原方程的形状参数
#YITA=np.exp(-b/k) #YITA为原方程的比例参数/强度参数
#b=-m*np.log(YITA)
return k*x+b
##偏差函数:x,y都是列表:这里的x,y跟上面的Xi,Yi中是一一对应的
def error(p,x,y):
return func(p,x)-y
p0=[1,20] #k,b的初始值
Para=leastsq(error,p0,args=(Xj,Yj))
#读取结果
k,b=Para[0]
YITA=np.exp(-b/k)
print("k=",k,"YITA=",YITA)
print("cost:"+str(Para[1]))
print("求解的拟合直线为:")
print("y="+str(round(k,2))+"x"+str(round(b,2)))
#画样本点
plt.figure(figsize=(8,6)) ##图像比例: 8:6
plt.scatter(Xj,Yj,color="green",label="样本数据",linewidth=2)
#画拟合直线
x=np.linspace(-10,10,100) ##在0-15直接画100个连续点
y=k*x+b ##函数式
plt.plot(x,y,color="red",label="拟合直线",linewidth=2)
plt.legend(loc='lower right') #绘制图例
plt.show()