最小二乘法的python的实现
偶然一次做项目的时候,需要对生长模型进行建模,这就意味着需要对其中的参数进行求解,其中最常见的也就是最小二乘法了,本文将其模型建立的过程的代码展示如下
第一步:导入关键的库
import pandas as pd
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
建立一个逻辑生长模型的类
1,首先初始化定义的类
def __init__(self, cases):
'''
Initializes Object
Args: 1-D array of cases at each time step
'''
self.parameters = np.random.exponential(size = 3)
print("Parameter Initialization")
print(self.parameters)
self.x = np.array([i for i in range(len(cases))])
self.y = np.array(cases)
2,定义需要拟合的方程
def logistic(self, t, a, b, c):
'''
Logistic function for training
Args: The time of the logistic prediction and the parameters
Returns: Output of logistic function
'''
return c / (1 + a * np.exp(-b*t))
3,curve_fit最小二乘法进行参数的求解
def trainLogistic(self):
'''
Trains logistic growth model
'''
bounds = (0, [1e10, 10, 1e10])
self.parameters, covariance = curve_fit(self.logistic, self.x, self.y, bounds=bounds, p0=self.parameters)
4,定义预测的函数,将求解后的方程带入进行求解
def predict(self, t):
'''
Logistic function for graphing and predictions
Args: The time of the logistic prediction
Returns: Output of logistic function
'''
return self.parameters[2] / (1 + self.parameters[0] * np.exp(-self.parameters[1]*t))
5,绘制求解的结果
def graph(self):
'''
Graphs the data with logistic model
'''
plt.scatter(self.x, self.y)
predictArr = np.vectorize(self.predict)
graphX = np.append(self.x, [i for i in range(len(self.x), len(self.x)*2)])
numOfDays = 0
for x in range(len(graphX[len(self.x):])):
if self.predict(graphX[len(self.x) + x]) > self.parameters[2]*(0.999):
numOfDays = x + 1
break
if (numOfDays == 0):
numOfDays = "More Than " + str(len(self.x))
plt.plot(graphX, predictArr(graphX))
plt.title('Logistic Model Predictions | Max at ' + str("%.1f" % self.parameters[2]) + "\nReached in " + str(numOfDays) + " Years")
plt.ylabel("Number")
plt.xlabel("Years")
print("Final Parameters")
print(self.parameters)
plt.show()
全部的代码如下所示:
class LogisticModel():
def __init__(self, cases):
'''
Initializes Object
Args: 1-D array of cases at each time step
'''
self.parameters = np.random.exponential(size = 3)
print("Parameter Initialization")
print(self.parameters)
self.x = np.array([i for i in range(len(cases))])
self.y = np.array(cases)
def logistic(self, t, a, b, c):
'''
Logistic function for training
Args: The time of the logistic prediction and the parameters
Returns: Output of logistic function
'''
return c / (1 + a * np.exp(-b*t))
def trainLogistic(self):
'''
Trains logistic growth model
'''
bounds = (0, [1e10, 10, 1e10])
self.parameters, covariance = curve_fit(self.logistic, self.x, self.y, bounds=bounds, p0=self.parameters)
def predict(self, t):
'''
Logistic function for graphing and predictions
Args: The time of the logistic prediction
Returns: Output of logistic function
'''
return self.parameters[2] / (1 + self.parameters[0] * np.exp(-self.parameters[1]*t))
def graph(self):
'''
Graphs the data with logistic model
'''
plt.scatter(self.x, self.y)
predictArr = np.vectorize(self.predict)
graphX = np.append(self.x, [i for i in range(len(self.x), len(self.x)*2)])
numOfDays = 0
for x in range(len(graphX[len(self.x):])):
if self.predict(graphX[len(self.x) + x]) > self.parameters[2]*(0.999):
numOfDays = x + 1
break
if (numOfDays == 0):
numOfDays = "More Than " + str(len(self.x))
plt.plot(graphX, predictArr(graphX))
plt.title('Logistic Model Predictions | Max at ' + str("%.1f" % self.parameters[2]) + "\nReached in " + str(numOfDays) + " Years")
plt.ylabel("Number")
plt.xlabel("Years")
print("Final Parameters")
print(self.parameters)
plt.show()