本示例中，x表示生产产品个数，Y表示生产成本。已知20天的生产数据，现要挖掘玩偶数量和生产成本关系。
一、导入数据，并可视化展出，根据图形判断选用模型。
image
代码如下：
import pandas as pd
import os
import matplotlib.pyplot as plt
def readData(dataPath):
"""
使用pandas读取数据
"""
data = pd.read_csv(dataPath)
return data
if __name__ == "__main__":
homePath = os.path.dirname(os.path.abspath(__file__))
# Windows下的存储路径与Linux并不相同
if os.name == "nt":
dataPath = "%s\\file\\test_goods.csv" % homePath
else:
dataPath = "%s/file/file/test_goods.csv" % homePath
test_goods_data = readData(dataPath)
plt.scatter(test_goods_data['x'],test_goods_data['y'])
#支持中文和负号
plt.rcParams["font.sans-serif"] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.title('散点图')
plt.xlabel('生产个数')
plt.ylabel('生产成本')
plt.show()
散点图如下：
image
根据散点图确定解决方向。目前由两个方向来确定。分别为机器学习和统计学。
二、 机器学习思路
2.1 确定问题类型
因为预测成本 是一个连续变量，而非离散变量，所以该问题为回归问题
我们搭建模型的目的，就是使模型尽可能精准的和真实 数据接近。
定义损失函数，使函数值最小 为预测值。
其实就是求 L 的最小值，至于为何用平方，是为了方便求导，取极值。
2.2 提取特征值
1、需要检查一下，数据中是否有异常，或者错误、为空等数据，如有进行删除或者补充等处理
2、对不能直接进行算数运算的数据，进行数学变换，转换成计算机课识别的数学语言。
3、对于可以直接进行数学运算的数据，也可以根据具体情况进行特征变换，比如平方，开根号，绝对值等等。
本例中不涉及特征转换。
2.3 确定模型形式，并估计参数。
由散点图可知，模型大致呈线性关系。可以直接使用线性模型。
我们的目的就是获取到参数 a,b的估计值 的值，使得 L 取得极小值。
2.4 决定系数
决定系数表明了，我们预测模型不能被解释的数据变化，占整体成本真实变化的的比例。
其中
2.5 代码如下
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn import linear_model
def evaluateModel(model, testData, features, labels):
"""
计算线性模型的均方差和决定系数
参数
----
model : LinearRegression, 训练完成的线性模型
testData : DataFrame，测试数据
features : list[str]，特征名列表
labels : list[str]，标签名列表
返回
----
error : np.float64，均方差
score : np.float64，决定系数
"""
# 均方差(The mean squared error)，均方差越小越好
error = np.mean(
(model.predict(testData[features]) - testData[labels]) ** 2)
# 决定系数(Coefficient of determination)，决定系数越接近1越好
score = model.score(testData[features], testData[labels])
return error, score
def visualizeModel(model, data, features, labels, error, score):
"""
模型可视化
"""
# 为在Matplotlib中显示中文，设置特殊字体
plt.rcParams['font.sans-serif']=['SimHei']
# 创建一个图形框
fig = plt.figure(figsize=(6, 6), dpi=80)
# 在图形框里只画一幅图
ax = fig.add_subplot(111)
# 在Matplotlib中显示中文，需要使用unicode
ax.set_title(u'%s' % "线性回归示例")
ax.set_xlabel('$x$')
# 转换Y轴注释方向
ax.set_ylabel('$y$',rotation='horizontal')
# 画点图，用蓝色圆点表示原始数据
ax.scatter(data[features], data[labels], color='b',
label=u'%s: $y = x + \epsilon$' % "真实值")
# 根据截距的正负，打印不同的标签
if model.intercept_ > 0:
# 画线图，用红色线条表示模型结果
ax.plot(data[features], model.predict(data[features]), color='r',
label=u'%s: $y = %.3fx$ + %.3f'\
% ("预测值", model.coef_, model.intercept_))
else:
ax.plot(data[features], model.predict(data[features]), color='r',
label=u'%s: $y = %.3fx$ - %.3f'\
% ("预测值", model.coef_, abs(model.intercept_)))
legend = plt.legend(shadow=True)
legend.get_frame().set_facecolor('#6F93AE')
# 显示均方差和决定系数
ax.text(0.99, 0.01,
u'%s%.3f\n%s%.3f'\
% ("均方差：", error, "决定系数：", score),
style='italic', verticalalignment='bottom', horizontalalignment='right',
transform=ax.transAxes, color='m', fontsize=13)
# 展示上面所画的图片。图片将阻断程序的运行，直至所有的图片被关闭
# 在Python shell里面，可以设置参数"block=False"，使阻断失效。
plt.show()
def trainModel(trainData, features, labels):
"""
利用训练数据，估计模型参数
参数
----
trainData : DataFrame，训练数据集，包含特征和标签
features : 特征名列表
labels : 标签名列表
返回
----
model : LinearRegression, 训练好的线性模型
"""
# 创建一个线性回归模型
model = linear_model.LinearRegression()
# 训练模型，估计模型参数
model.fit(trainData[features], trainData[labels])
return model
def linearModel(data):
"""
线性回归模型建模步骤展示
参数
----
data : DataFrame，建模数据
"""
features = ["x"]
labels = ["y"]
# 划分训练集和测试集
trainData = data[:15]
testData = data[15:]
# 产生并训练模型
model = trainModel(trainData, features, labels)
# 评价模型效果
error, score = evaluateModel(model, testData, features, labels)
# 图形化模型结果
visualizeModel(model, data, features, labels, error, score)
def readData(path):
"""
使用pandas读取数据
"""
data = pd.read_csv(path)
return data
if __name__ == "__main__":
homePath = os.path.dirname(os.path.abspath(__file__))
# Windows下的存储路径与Linux并不相同
if os.name == "nt":
dataPath = "%s\\file\\test_goods.csv" % homePath
else:
dataPath = "%s/file/test_goods.csv" % homePath
data = readData(dataPath)
linearModel(data)
训练结果如图：
image
2.6 总结
具体操作步骤如下：
1、将数据集划分为测试集和训练集
2、利用训练集，训练线性模型，估计相关参数。
3、利用测试集，评估模型效果，这里，主要用了两个参数度。均方差和决定系数。
这部分的实现过程，过于机械。而且也没有优化步骤。只是机械的根据数据集来训练。
三、统计方法代码实现
3.1 首先我们已经由图形化的数据，看出是线性回归模型，那么就可以利用 Statsmodels 进行训练。
statsmodels 这个Python库是集合了多种统计模型，可执行统计测试以及数据探索可视化。主要包含的子模块有：
1、回归模型：线性回归 ，通用线性回归，鲁邦线性模型 ，线性混合效应模型等。
2、方差分析（ANOVA）。
3、时间序列分析：AR , ARMA , ARIMA , VAR等。
4、非参数方法： 核密度估计 ， 核回归。
5、统计模型结果可视化。
3.2 代码如下
import os
import sys
import numpy as np
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std
import matplotlib.pyplot as plt
import pandas as pd
def modelSummary(re):
"""
分析线性回归模型的统计性质
"""
# 整体统计分析结果
print(re.summary())
# 在Windows下运行此脚本需确保Windows下的命令提示符(cmd)能显示中文
# 用f test检测x对应的系数a是否显著
print("检验假设x的系数等于0：")
print(re.f_test("x=0"))
# 用f test检测常量b是否显著
print("检测假设const的系数等于0：")
print(re.f_test("const=0"))
# 用f test检测a=1, b=0同时成立的显著性
print("检测假设x的系数等于1和const的系数等于0同时成立：")
print(re.f_test(["x=1", "const=0"]))
def visualizeModel(re, data, features, labels):
"""
模型可视化
"""
# 计算预测结果的标准差，预测下界，预测上界
prstd, preLow, preUp = wls_prediction_std(re, alpha=0.05)
# 为在Matplotlib中显示中文，设置特殊字体
plt.rcParams['font.sans-serif']=['SimHei']
# 创建一个图形框
fig = plt.figure(figsize=(6, 6), dpi=80)
# 在图形框里只画一幅图
ax = fig.add_subplot(111)
# 在Matplotlib中显示中文，需要使用unicode
ax.set_title(u'%s' % "线性回归统计分析示例")
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
# 画点图，用蓝色圆点表示原始数据
ax.scatter(data[features], data[labels], color='b',label=u'%s: $y = x + \epsilon$'%"真实值")
# 画线图，用红色虚线表示95%置信区间
ax.plot(data[features], preUp, "r--", label=u'%s' % "95%置信区间")
ax.plot(data[features], re.predict(data[features]), color='r',
label=u'%s: $y = %.3fx$'\
% ("预测值", re.params[features]))
ax.plot(data[features], preLow, "r--")
legend = plt.legend(shadow=True)
legend.get_frame().set_facecolor('#6F93AE')
plt.show()
def trainModel(X, Y):
"""
训练模型
"""
model = sm.OLS(Y, X)
re = model.fit()
return re
def linearModel(data):
"""
线性回归统计性质分析步骤展示
参数
----
data : DataFrame，建模数据
"""
features = ["x"]
labels = ["y"]
Y = data[labels]
# 加入常量变量
X = sm.add_constant(data[features])
# 构建模型
re = trainModel(X, Y)
# 分析模型效果
modelSummary(re)
# const并不显著，去掉这个常量变量
resNew = trainModel(data[features], Y)
# 输出新模型的分析结果
print(resNew.summary())
# 将模型结果可视化
visualizeModel(resNew, data, features, labels)
def readData(path):
"""
使用pandas读取数据
"""
data = pd.read_csv(path)
return data
if __name__ == "__main__":
homePath = os.path.dirname(os.path.abspath(__file__))
# Windows下的存储路径与Linux并不相同
if os.name == "nt":
dataPath = "%s\\file\\test_goods.csv" % homePath
else:
dataPath = "%s/file/test_goods.csv" % homePath
data = readData(dataPath)
linearModel(data)
3.3 输出如下：
image
上方参数分析如下：
- Dep.Variable:输出变量名称
- Model: 模型名称
- Method: 方法 （least squares就表示最小二乘法）
- Date,Time 表示时间日期
- No.Observations: 观察的数据量（样本数据量）
- Df Residuals: 残差自由度
- Df model : 模型参数个数，也就是自变量个数
-
R- squared : 可决系数（决定系数），用来判断估计的准确性，范围在 [0,1] 约接近1 ，说明对y的解释能力越强，拟合越好。解释如下。
image
- Adj-R- squared: 通过样本数量与模型数量对R-squared进行修正，奥卡姆剃刀原理，避免描述冗杂
- F-statistic : 衡量拟合的显著性, 重要程度。模型的均方误差除以残差的均方误差,值越大,越不可能
- Prob(F-statistic): 当prob（F-statistic）<α时，表示拒绝原假设，即认为模型是显著的； 当prob（F-statistic）>α时，表示接受原假设，即认为模型不是显著的
- Log likelihood （对数似然比LLR） :(很多说法是值越大，说明模型拟合的较好)(但有待考察,似然比服从统计量,大于卡方临界值拒绝原假设)
- AIC: AIC可以表示为： AIC=2k-2ln(L) 其中：k是参数的数量，L是似然函数。 衡量拟合优良性,选择AIC 最小的模型, 引入了惩罚项,避免参数过多,过拟合
- BIC: 贝叶斯信息准则 BIC=kln(n)-2ln(L) ,BIC相比AIC在大数据量时对模型参数惩罚得更多，导致BIC更倾向于选择参数少的简单模型。
- coef: 系数 const表示常数项
- std err :系数估计的基本标准误差
- t : t 统计值,衡量系数统计显著程度的指标
- P>|t| : 系数= 0的零假设为真的P值。如果它小于置信水平，通常为0.05，则表明该术语与响应之间存在统计上显着的关系。度的指标
- 区间[0.025,0.975]: 95％置信区间的下限和上限值
- Omnibus :属于一种统计测验,测试一组数据中已解释方差是否显著大于未解释方差,但omnibus不显著,模型也可能存在合法的显著影响, 比如两个变量中有一个不显著,即便另一个显著.通常用于对比
- Prob(Omnibus):将上面的统计数据变成概率
- Durbin-Watson : 残差是否符合正态分布,在2左右说明是服从正态分布的,偏离2太远,解释能力受影响
是否自相关, 受到前后影响 ,与表中上限进行比较,如果D>上限 不存在相关性 .
D<下限 存在正相关性,在上下限之间,无法得出结论 - Skewness: 偏度, 关于平均值的数据对称性的度量。正态分布误差应是关于平均值对称的分布。
- Kurtosis: 峰度, 分布形状的量度,比较接近均值与远离均值的数据量 如果大于三,说明峰的形状比较陡峭,形状较尖
正态分布的峰度（系数）为常数3，均匀分布的峰度（系数）为常数1.8 - Jarque-Bera(JB) : Jarque–Bera检验是对样本数据是否具有符合正态分布的偏度和峰度的拟合优度的检验。
其统计测试结果总是非负的。如果结果远大于零，则表示数据不具有正态分布。 - Prob(JB): 上面统计量的概率形式
- Cond. No :多重共线性测试(如果多个参数,这些参数是否相互关联)
3.4 输出分析
由输出结果可知，b可以舍去，因为其等于0 时的概率为32.3% ，并不显著。
所以模型改写后为：
网友评论