玩具生产成本预测模型-优快云博客

本文链接：https://blog.youkuaiyun.com/techfei/article/details/99704268

问题背景，玩具生产个数与生产成本的关系如下：

生产记录
日期	玩具数量	成本	天数
04/01	10	7.7	1
04/02	10	9.87	2
04/03	11	10.87	3
04/04	12	12.18	4
04/05	13	11.43	5
04/06	14	13.36	6
04/07	15	15.15	7
04/08	16	16.73	8
04/09	17	17.4	9
...	...	...	...

分析玩具数量与成本之间的关系。

1 从机器学习角度看线性回归

确定场景类型

（1）在数据里有需要被预测的变量：成本，所以这是一个监督式学习。

（2）被预测的变量是一个连续的量，所以这是一个回归问题。

定义损失函数

（1）搭建模型的目标是让预测值和实际成本接近

（2）损失函数为

$LL=\sum _{i}\left | y_{i}- \hat{y_{i}}\right |$

由于这个公式不可导，所以换成一个处处可导的函数如下，

特征提取

（1）去掉错误和异常数据

（2）变量本身的数学运算（加减乘除）有意义，就可以在模型中使用

（3）对变量做数学变换，比如平方，可以得到新的特征。

确定模型并估计参数

选用线性模型

$\hat{y_{i}}=ax_{i}+b$

a表示生产一个玩具的变动成本，b表示固定成本。

参数的估计值 $(\hat{a},\hat{b})$ 应该使损失函数达到最小值

$(\hat{a},\hat{b})=argmin_{a,b}\sum _{i}(y_{i}-ax_{i}-b)^{2}$

评估模型效果

（1）定义线性模型的均方差

$MSE=\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\hat{y_{i}})^{2}=\frac{1}{n}L$

均方差越小，效果越好。

（2）从数据解释性上说，我们希望未被模型解释的成本 $(y_{i}-\hat{y_{i}})$ 占成本变化 $(y^{_{i}}-\bar{y})$ 的比例越小越好，由此定义决定系数，决定系数越接近1越好

$\bar{y}=\frac{1}{n}\sum y_{i}$

$SS_{tot}=\sum (y_{i}-\bar{y})^{2}$

$SS_{res}=\sum (y_{i}-\hat{y_{i}})^{2}$

$R^{2}=1-\frac{SS_{res}}{SS_{tot}}$

2 从统计学角度看线性回归

假设条件概率

假设生产成本 $y_{i}$ 和生产数量 $x_{i}$ 之间的关系如下

$y_{i}=ax_{i}+b+\varepsilon _{i}$

其中 $\varepsilon_{i} \sim N(0,\sigma ^{2})$ ,当a,b,x都确定时， $y_{i}\sim N(ax_{i}+b,\sigma ^{2})$ , $\begin{Bmatrix} y_{i} \end{Bmatrix}$ 相互独立，即条件概率满足

$P(y_{i}|a,b,x_{i},\sigma ^{2})\sim N(ax_{i}+b,\sigma ^{2})$

估计参数

$\begin{Bmatrix} y_{i} \end{Bmatrix}$ 出现的联合概率（似然函数）如下

$\ln P(Y|a,b,X,\sigma ^{2})=-0.5n\ln (2\pi \sigma ^{2})-(1/2\sigma ^{2})\sum (y_{i}-ax_{i}-b)^{2}$

使这个概率最大的参数是参数估计的最佳选择，因此这种方法叫最大似然估计（MLE）。

参数（a,b）的估计值如下

$(\hat{a},\hat{b})=argmax_{a,b}\ln P(Y|a,b,X,\sigma ^{2})=argmin_{a,b}\sum (y_{i}-ax_{i}-b)^{2}$

同理可得参数 $\sigma ^{2}$ 估计值如下

$\hat{\sigma ^{2}}=argmax_{a,b}P(Y|a,b,X,\sigma ^{2})=\sum (y_{i}-\hat{y_{i}})^{2}/n$

推导参数分布

参数估计值 $(\hat{a},\hat{b},\hat{\sigma ^{2}})$ 都是随机变量，服从均值为参数真实值的正态分布

$\hat{a}\sim N(a,\sigma ^{2}/\sum (x_{i}-\bar{x})^{2})$

$\hat{b}\sim N(b,\sigma ^{2}/n)$

$\hat{\sigma ^{2}}\sim X_{n-2\sigma ^{2}/n}^{2}$

参数估计值的方差随着数据量的增大而减少，即数据量越大，模型估计的参数就越接近真实值。

假设检验与置信区间

置信区间：对于被预测对象，真实值的大致范围是怎样

P（c1<=u<=c2）=1-a

其中，a是显著性水平，1-a是置信水平

假设检验：参数等于0的概率是否小于1%

3 模型实现

机器学习代码实现

homePath = os.path.dirname(os.path.abspath(__file__))

os.path.abspath(__file__)：获取当前脚本的完整路径，

homePath = os.path.dirname：去掉文件名，返回目录路径

用机器学习方法解决线性回归问题的代码实现里主要涉及四个函数：读取数据，训练模型，评价效果，可视化

# -*- coding: UTF-8 -*-
"""
此脚本用于展示使用sklearn搭建线性回归模型
"""


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
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.set_title(u'%s' % "线性回归示例")
    else:
        ax.set_title(u'%s' % "线性回归示例".decode("utf-8"))
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    # 画点图，用蓝色圆点表示原始数据
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.scatter(data[features], data[labels], color='b',
            label=u'%s: $y = x + \epsilon$' % "真实值")
    else:
        ax.scatter(data[features], data[labels], color='b',
            label=u'%s: $y = x + \epsilon$' % "真实值".decode("utf-8"))
    # 根据截距的正负，打印不同的标签
    if model.intercept_ > 0:
        # 画线图，用红色线条表示模型结果
        # 在Python3中，str不需要decode
        if sys.version_info[0] == 3:
            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'\
                % ("预测值".decode("utf-8"), model.coef_, model.intercept_))
    else:
        # 在Python3中，str不需要decode
        if sys.version_info[0] == 3:
            ax.plot(data[features], model.predict(data[features]), color='r',
                label=u'%s: $y = %.3fx$ - %.3f'\
                % ("预测值", model.coef_, abs(model.intercept_)))
        else:
            ax.plot(data[features], model.predict(data[features]), color='r',
                label=u'%s: $y = %.3fx$ - %.3f'\
                % ("预测值".decode("utf-8"), model.coef_, abs(model.intercept_)))
    legend = plt.legend(shadow=True)
    legend.get_frame().set_facecolor('#6F93AE')
    # 显示均方差和决定系数
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        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)
    else:
         ax.text(0.99, 0.01, 
            u'%s%.3f\n%s%.3f'\
            % ("均方差：".decode("utf-8"), error, "决定系数：".decode("utf-8"), 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\\data\\simple_example.csv" % homePath
    else:
        dataPath = "%s/data/simple_example.csv" % homePath
    data = readData(dataPath)
    linearModel(data)

统计学代码实现

t-检验只能针对一个随机变量，f-检验可以一个或多个随机变量

搭建模型是一个需要根据数据和模型结果，不断修正，迅速迭代的过程

# -*- coding: UTF-8 -*-
"""
此脚本用于展示如何使用statsmodels搭建线性回归模型
"""


# 保证脚本与Python3兼容
from __future__ import print_function

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
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.set_title(u'%s' % "线性回归统计分析示例")
    else:
        ax.set_title(u'%s' % "线性回归统计分析示例".decode("utf-8"))
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    # 画点图，用蓝色圆点表示原始数据
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.scatter(data[features], data[labels], color='b', 
            label=u'%s: $y = x + \epsilon$' % "真实值")
    else:
        ax.scatter(data[features], data[labels], color='b',
            label=u'%s: $y = x + \epsilon$' % "真实值".decode("utf-8"))
    # 画线图，用红色虚线表示95%置信区间
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        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]))
    else:
        ax.plot(data[features], preUp, "r--", label=u'%s' % "95%置信区间".decode("utf-8"))
        ax.plot(data[features], re.predict(data[features]), color='r',
            label=u'%s: $y = %.3fx$'\
            % ("预测值".decode("utf-8"), 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\\data\\simple_example.csv" % homePath
    else:
        dataPath = "%s/data/simple_example.csv" % homePath
    data = readData(dataPath)
    linearModel(data)