7. 统计回归模型

数学建模

张瑞
中国科学技术大学数学科学学院

统计回归模型

此类模型的建模思想是:

  • 模型通常涉及大量的数据;实际问题是要从这些数据中提取有用信息。
  • 需要利用统计学知识,对数据进行量化分析、参数估计和假设检验等处理,从中得到相对可信赖的 信息。
  • 回归分析是研究多个因素之间的定量关系的最常用的统计方法。

数理统计的基本知识

数学建模过程中，会遇到：有一项需要研究的指标$X$（如：产品质量，农作物产量等）受到一些因素$A_1$, $A_2$, $\cdots$, $A_m$的影响（如：设备、原材料、工艺对产品质量的影响，种子、肥料、农药对农作物产量的影响），此外$X$还受到一些随机误差的影响。我们想了解因素$A_1$, $A_2$, $\cdots$, $A_m$对指标$X$的影响是否显著，随机误差对$X$的影响有多大。

  • 数理统计(简称统计)研究的对象是受随机因素影响的数据,是以概率论为基础的一门应用学科。
  • 面对一批数据如何进行描述与分析,需要掌握参数估计和假设检验这两个数理统计的最基本方法。

统计上要研究对象的全体为总体。总体中的一个随机样品称为个体。从总体中随机产生的若干个个体的集合称为样本。样本中个体的数量称为样本容量。为了获取样本，需要对总体中的个体进行抽样试验，简称抽样。

简单地说,统计的任务是由样本推断总体。具体来说，根据某种原则或经验假设总体服从于某种概率分布，而通过样本推断总体所服从的概率分布的若干参数。

设$X$为一个总体，分布函数为$F(x)$，$x=(x_1$,$\cdots$,$x_n$为从总体中取出的容量为$n$的样本。

  • 这引些样本数据需要对它进行一定的加工,才能提出有用的信息,用作对总体(分布)参数的估计和检验。

统计量就是加工出来的、反映样本数量特征的函数,它不含任何未知量。

  • 表示位置的统计量(算术平均值和中位数)
    1. 算术平均值(简称均值)描述数据取值的平均位置,记作 $\bar x$ ,
      \[\bar x=\frac1n\sum_{i=1}^nx_i \]
    2. 中位数是将数据由小到大排序后位于中间位置的那个数值。

Matlab 

mean(x)返回 x 的均值,median(x)返回中位数。

Python 

import statistics as st

mean() 数据的算数平均值
median() 数据的中位数

numpy.mean()
numpy.median()
numpy.ptp()
scipy.stats.scoreatpercentile(a, per[, limit, ...])

Mathematica 

(* 用 Mean 求 data 的均值： *)
Mean[data]
(* 用 Median 求中位数 *)
Median[data]
  • 表示变异程度的统计量(标准差、方差和极差)
    1. 标准差 $s$ 定义为
      \[s=\left[\frac1{n-1}\sum_{i=1}^n(x_i-\bar x)^2\right]^{\frac12} \]
      它是各个数据与均值偏离程度的度量,这种偏离不妨称为变异
    2. 方差是标准差的平方$s^2$
    3. 极差是 $x = ( x_1 , \cdots , x_n )$ 的最大值与最小值之差。

Matlab 

std(x)返回 x 的标准差,var(x)返回方差,range(x)返回极差。

Python 

import statistics as st

stdev() 数据样本的标准差
variance() 数据样本的方差

numpy.std()
numpy.var()

Mathematica 

(* 可以求 Variance 和 StandardDeviation： *)
Variance[data]
StandardDeviation[data]
(* 如果您想弄清楚分位数，可以使用 Quantile. 第一个参数是一个数据集合；第二个参数是一个0到1之间的数 q：*)
Quantile[data, 4/5]
  • 表示分布形状的统计量(偏度和峰度)
    1. 随机变量 $x$ 的偏度 $\displaystyle v_1=\frac1{s^3}\sum_{i=1}^n(x_i-\bar x)^3$，它反映分布的对称性
    2. 峰度指的是 $\displaystyle v_2=\frac1{s^4}\sum_{i=1}^n(x_i-\bar x)^4$，可以用作衡量偏离正态分布的尺度之一(正态分布的峰度为3，峰度大，说明样本中含有较多远离均值的数据)。
    3. $k$阶中心距 $D_k=\frac1n\sum_{i=1}^n(x_i-\bar x)^k$, $k=1,2,\cdots$

Matlab 中 

moment(x,order)返回 x 的 order 阶中心矩,order 为中心矩的阶数。
skewness(x)返回 x 的偏度,
kurtosis(x)返回峰度。

Python 

scipy.stats.skew
scipy.stats.kurtosis (返回0表示是正态分布)
scipy.stats.norm (中心距)

Mathematica 

(* 给出 list 中元素的峰度系数. *)
Kurtosis[list]
(* 给出列表 list 中元素的偏度系数. *)
Skewness[list]
(* 给出 list 中元素相应均值的 r 阶中心矩. *)
CentralMoment[list,r]

随机变量 $x$ 的 $r$ 阶中心矩为 $E( x − Ex )$ 。

\[\begin{aligned} v_1=E\left[\left(\frac{x-E(x)}{\sqrt{Dx}}\right)^3\right]=\frac{E\left[(x-Ex)^3\right]}{(Dx)^{\frac{3}{2}}} \\ v_2=E\left[\left(\frac{x-E(x)}{\sqrt{Dx}}\right)^4\right]=\frac{E\left[(x-Ex)^4\right]}{(Dx)^{2}} \\ \end{aligned} \]
  • 统计中几个重要的概率分布

正态分布 正态分布随机变量 $X$ 的密度函数曲线呈中间高两边低、对称的钟形,期望(均值)$EX = \mu$ ,方差 $DX = \sigma^2$ ,记作 $X \sim N (\mu , \sigma^2)$ , $\sigma$ 称均方差或标准差,当 $\mu = 0$ , $\sigma = 1$时称为标准正态分布,记作 $X \sim N ( 0 , 1 )$ 。正态分布完全由均值 $\mu$ 和方差 $\sigma$ 决定,它的偏度为 0,峰度为 3。

  • 正态分布可以说是最常见的(连续型)概率分布
  • 中心极限定理表明，在大量相互独立的、作用差不多大的随机因素影响下形成的随机变量,其极限分布为正态分布
  • 3 个数字
    • 68%的数值落在距均值左右 1 个标准差的范围内,即
      \[P ( \mu-\sigma\leq X\leq\mu+\sigma )=0.68 \]
    • 95%的数值落在距均值左右 2 个标准差的范围内,即
      \[P (\mu − 2 \sigma \leq X \leq \mu+2\sigma)=0.95 \]
    • 99.7%的数值落在距均值左右 3 个标准差的范围内,即
      \[P (\mu − 3 \sigma \leq X\leq \mu+3\sigma)=0.97 \]

  • $\chi^2$ 分布(Chi square) 若$X_1$,$\cdots$,$X_n$相互独立、服务标准正态分布$N(0,1)$的随机变量，则它们的平方和$Y =\sum_{i=1}^nX_i^2$服从 $\chi^2$ 分布,记作 $Y\sim \chi^2(n)$，$n$为自由度。它的期望$EY=n$，方差$DY=2n$。

  • $t$分布 若$X\sim N(0,1)$，$Y\sim \chi^2(n)$，且相互独立，则$T=\frac{X}{\sqrt{Y/n}}$服从$t$分布，记为$T\sim t(n)$，$n$为自由度。$t$分布又称学生氏(Student)分布。

  • $F$分布 若$X\sim \chi^2(n_1)$，$Y\sim\chi^2(n_2)$，且相互独立，则$F=\frac{X/n_1}{Y/n_2}$服从$F$分布，记为$F\sim F(n_1,n_2)$，$(n_1, n_2)$称自由度。

Matlab 

norm- 正态分布， chi2- Chi Square
t - t分布， f- F分布

函数 

pdf 概率密度; cdf 分布函数; inv 分布函数的反函数;
stat 均值与方差; rnd 随机数生成

如： 

p=normpdf(x,mu,sigma) 均值 mu、标准差 sigma 的正态分布在 x 的密度函数
                      (mu=0,sigma=1 时可缺省)。
p=tcdf(x,n) t 分布(自由度 n)在 x 的分布函数。
x=chi2inv(p,n) Chi Square 分布(自由度 n)使分布函数 F(x)=p 的 x(即 p 分位数)。
[m,v]=fstat(n1,n2) F 分布(自由度 n1,n2)的均值 m 和方差 v。

Mathematica 

(* 正态分布的 概率密度与分布 函数 *)
PDF[NormalDistribution[mu,sigma],x]
CDF[NormalDistribution[mu,sigma],x]
(* 一些分布 *)
ChiSquareDistribution[mu]
FRatioDistribution[n,m]
StudentTDistribution[v]
(* https://reference.wolfram.com/language/tutorial/ContinuousDistributions.html *)

Python 

scipy.stats.norm
scipy.stats.t
scipy.stats.lognormal
# 可以接入方法
rvs(..) # 产生相应分布的随机分布数列
pdf(..) # 给出概率密度函数
# 如
scipy.stats.norm.pdf(..)
scipy.stats.t.pdf(..)
  • 用样本来推断总体,需要知道样本统计量的分布,而样本又是一组与总体同分布的随机变量,所以样本统计量的分布依赖于总体的分布。
  • 当总体服从一般的分布时,求某个样本统计量的分布是很困难的,只有在总体服从正态分布时,一些重要的样本统计量(均值、标准差)的分布才有便于使用的结果。
  • 另一方面,现实生活中需要进行统计推断的总体,多数可以认为服从(或近似服从)正态分布,所以统计中人们在正态总体的假定下研究统计量的分布,是必要的与合理的。

设总体 $X \sim N ( \mu, \sigma)$ , $x_1$ , $\cdots$, $x_n$ 为一容量 $n$ 的样本,则用其均值 $\bar x$ 和标准差 $s$ 构造的下面几个分布在统计中是非常有用的。

\[\bar x \sim N(\mu, \frac{\sigma^2}{n}) ,\mbox{or} \ \ \frac{\bar x-\mu}{\sigma/\sqrt{n}}\sim N(0,1) \]
\[\frac{(n-1)s^2}{\sigma^2}\sim\chi^2(n-1) \]
\[\frac{\bar x-\mu}{s/\sqrt{n}}\sim t(n-1) \]

设有$X\sim N(\mu_1,\sigma_1^2)$和$Y\sim N(\mu_2,\sigma_2^2)$，和$\bar x$, $\bar y$, $s_1$, $s_2$，则

\[\frac{(\bar x-\bar y)-(\mu_1-\mu_2)}{\sqrt{\sigma_1^2/n_1+\sigma_2^2/n_2}}\sim N(0,1) \]
\[\frac{(\bar x-\bar y)-(\mu_1-\mu_2)}{s_w\sqrt{1/n_1+1/n_2}}\sim t(n_1+n_2-2) \]

其中$s_w^2=\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}$

参数估计

利用样本对总体进行统计推断的一类问题是参数估计,即假定已知总体的分布,通常是 $X \sim N (\mu, \sigma)$,估计有关的参数,如 $\mu$, $\sigma$。参数估计分点估计和区间估计两种。

  • 点估计 是用样本统计量确定总体参数的一个数值。评价估计优劣的标准有无偏性、最小方差性、有效性等,估计的方法有矩法、极大似然法等。

最常用的是对总体均值 $\mu$ 和方差 $\sigma^2$ (或标准差 $\sigma$ )作点估计。让我们暂时抛开评价标准,当从一个样本算出样本均值 $\bar x$ 和方差 $s$ 后,一个自然、合理的点估计显然是

\[\hat \mu=\bar x, \hat \sigma^2=s^2, \hat \sigma=s \]

(在字母上加^表示它的估计值)

Matlab 

[mu,sigma,muci,sigmaci]=normfit(x,alpha)

Python 

scipy.stats.norm.fit(data,loc=0,scale=1)
  • 同一个总体的2个样本，得到的均值往往是不同的。点估计没有告诉我们这个估计值的精度和可信程度。
  • 解决这类问题的办法是：以一定的置信概率确定待估参数（如总体均值）的一个区间，称为区间估计。

一般地,总体的待估参数记作 $\theta$ (如 $\mu$ 或 $\sigma^2$ ),由样本算出的 $\theta$ 的估计量记作 $\hat \theta$，人们常希望给出一个区间 $[ \hat \theta_1, \hat \theta_2]$,使 $\theta$ 以一定的概率落在此区间内。若有

\[P(\hat\theta_1<\theta<\hat\theta_2)=1-\alpha, 0<\alpha<1 \]

则 $[ \hat \theta_1, \hat \theta_2]$ 称为 $\theta$ 的置信区间, $\hat\theta_1$, $\hat\theta_2$分别称为置信下限和置信上限, $1-\alpha$ 称为置信概率或置信水平,$\alpha$ 称为显著性水平。

给出的置信水平为 $1-\alpha$ 的置信区间 $[ \hat \theta_1, \hat \theta_2]$ ,称为 $\theta$ 的区间估计。

  • 置信区间越小,估计的精度越高;置信水平越大,估计的可信程度越高。
  • 但是这两个指标显然是矛盾的,通常是在一定的置信水平下使置信区间尽量小。
  • 通俗地说,区间估计给出了点估计的误差范围。
  • 总体正态$N(\mu,\sigma^2)$的假定下，对总体均值$\mu$的区间估计
    1. 总体方差$\sigma^2$已知。（Z检验，或$u$估计）
    区间估计的基本思想是用样本构造一个已知其分布的统计量。设样本容量$n$，均值$\bar x$。则对正态总体$N(\mu,\sigma^2)$有$E\bar x=\mu$, $D\bar x=\sigma^2/n$，$\bar x\sim N(\mu,\sigma^2/n)$，则
    \[u=\frac{\bar x-\mu}{\sigma\sqrt n} \sim N(0,1) \]
    也就是说，当总体是正态分布$N(\mu,\sigma$时，$u$就应该满足$N(0,1)$。给定一个较小的数$\alpha$(取$0.05$或$0.01$)，取$N(0,1)$的$1-\alpha/2$分位数为$u_{1-\frac{\alpha}2}$，则应该有
    \[-u_{1-\frac{\alpha}2}\leq u\leq u_{1-\frac{\alpha}2} \]

$u$应该满足

\[P\left(-u_{1-\frac{\alpha}2}\leq u\leq u_{1-\frac{\alpha}2}\right)=1-\alpha \]

作变换后，有

\[P\left(\bar x-u_{1-\frac{\alpha}2}\frac{\sigma}{\sqrt n}\leq \mu\leq \bar x+u_{1-\frac{\alpha}2}\frac{\sigma}{\sqrt n}\right)=1-\alpha \]

即在置信水平$1-\alpha$下，$\mu$的置信区间为

\[\left[\bar x-u_{1-\frac{\alpha}2}\frac{\sigma}{\sqrt n}\leq \mu, x+u_{1-\frac{\alpha}2}\frac{\sigma}{\sqrt n}\right] \]

因此，粗略地来说：$1-\alpha=0.95$，重复取100个样本，得到的100个置信区间中约有95个包含了总体均值$\mu$。而我们只有1个样本，这个样本的置信区间中不含$\mu$的概率为$\alpha=0.05$。

随机变量的特性完全由它的分布函数($F(x)=P(\xi\leq x)$)或密度函数(满足$F(x)=\int_{-\infty}^xp(x)dx$的$p(x)$)来描述。对于$0<\alpha<1$，若$F(x_\alpha)=1-\alpha$，则称$x_\alpha$为这个分布的上$\alpha$分位数。

Python 

scipy.stat.norm.ppf(1-0.05) # 分位点\alpha=0.05

  1. 总体方差$\sigma^2$不知道时。

    只能用样本方差$s^2$来代替$\sigma^2$。但是$\frac{\bar x-\mu}{s/\sqrt n}$不再服从$N(0,1)$。可以证明
    \[t=\frac{\bar x-\mu}{s/\sqrt n}\sim t(n-1) \]
    $t$分布也是对称的。类似地，$t(n-1)$的$1-\alpha/2$分位数为$t_{1-\alpha/2}$，则在置信水平$1-\alpha$下，$\mu$的置信区间为
    \[\left[\bar x-t_{1-\alpha/2}\frac{s}{\sqrt n}, \bar x+t_{1-\alpha/2}\frac{s}{\sqrt n}\right] \]
  • 两个总体的均值 记两个独立的总体为$X\sim N(\mu_1, \sigma_1^2)$，$Y\sim N(\mu_2, \sigma_2^2)$。两个分别来自$X$和$Y$的样本的容量为$n_1$和$n_2$，样本均值为$\bar x$和$\bar y$，样本标准差$s_1$和$s_2$。
    1. 当$\sigma_1$, $\sigma_2$已知时，可以证明
      \[u=\frac{\bar x-\bar y}{\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}}\sim N(0,1) \]
    2. 当$\sigma_1$, $\sigma_2$未知时，有
      \[t=\frac{\bar x-\bar y}{\sqrt{\frac{s^2}{n_1}+\frac{s^2}{n_2}}}\sim t(n_1+n_2-2),\ \ s^2=\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2} \]

假设检验

统计推断的另一类重要问题是假设检验问题。

例 1. 甲方向乙方供货，约定废品率不超过3%，预先约定显著性水平$\alpha=0.05$。现在从一批产品中抽取100件，发现有5件废品，问乙方是否应该接受这批产品？

假设：这批产品是合格的，即至多3%的废品率。

检验：从总的产品中，取出100件产品，含5件废品的概率有多大？

  • 如果这是小概率事件，则假设不成立；
  • 否则，就可以认为假设是成立的。

若货物总量是$n\times 100$件，货物最多$3n$件废品。则取出100件中，含5件废品的概率为

\[P(n)=\frac{C^5_{3n}C_{97n}^{95}}{C_{100n}^{100}} \]

可以得到

\[\begin{aligned} P(2)=0.091, &P(30)=0.1017, \\ P(500)=0.10133, &P(8000)=0.10131 \end{aligned} \]

假设检验就是先对总体作出假设，再根据样本对假设做出判断：是接受还是拒绝。

  • 假定总体服务正态分布，称为参数假设。常用的检验方法有$z$检验法，$t$检验法，$F$检验法和$\chi^2$检验法
  • 假定总体服务其它分布（泊松分布）或某个未知分布，称为非参数检验。

参数检验(基于确认的总体分布)

  • 单个总体 $N ( \mu,\sigma)$ 均值 $\mu$ 的检验
    1. $\sigma^2$已知，关于$\mu$的检验（$Z$检验）
    2. $\sigma^2$未知，关于$\mu$的检验（$t$检验）
  • 两个正态总体均值差的检验( t 检验)

设$X\sim N(\mu,\sigma^2)$，$x_1$, $\cdots$, $x_n$为一个样本，且

\[S=\sqrt{\frac1{n-1}\sum_{i=1}^n(x_i-\bar x)^2}, S_0^2=\frac{1}n\sum_{i=1}^n(x_i-\mu)^2 \]

关于$\mu$及$\sigma^2$的假设检验为

假设$H_0$ 条件 统计量 统计量分布 不定域
$\mu=\mu_0$ 已知 $\sigma^2$ $z=\frac{\bar x-\mu_0}{\sigma}\sqrt{n}$ $N(0,1)$ $\|z\|>\alpha$
$\mu=\mu_0$ 未知 $\sigma^2$ $t=\frac{\bar x-\mu_0}{S}\sqrt{n}$ $t(n-1)$ $\|t\|>t_\alpha(n-1)$
$\sigma=\sigma_0^2$ 已知$\mu$ $\chi^2=\frac{nS_0^2}{\sigma_0^2}$ $\chi^2(n)$
$\sigma=\sigma_0^2$ 未知$\mu$ $\chi^2=\frac{(n-1)S^2}{\sigma_0^2}$ $\chi^2(n-1)$

非参数检验(总体不服从正态分布, 不对分布的参数来进行推断)

  • $\chi^2$检验（皮尔逊）
  • Wilcoxon秩和检验
  • 分布拟合检验: 在实际问题中,有时不能预知总体服从什么类型的分布,这时就需要根据样本来检验关于分布的假设。
    1. 偏度、峰度检验（Jarque-Bera检验，适用于大样本）
    2. Kolmogorov-Smirnov检验（与给定分布函数做比较，通常做标准正态性检验）
    3. Lilliefors检验（KS检验的改进用于一般正态性检验）
  • 中位数检验
  • 0-1分布总体均值的假设检验

如：一件产品要么是合格品，要么是废品（0或1）。总体服务0-1分布。若废品率为$p$，则$X$的期望为$\mu=p$，方差$\sigma^2=p(1-p)$，不服从正态分布。但当样本容量$n$足够大时，对样本均值$\bar x$，有$u=\frac{\bar x-\mu}{\sigma/\sqrt n}$近似地服从$N(0,1)$，即

\[u=\frac{\bar x-p_0}{\sqrt{p_0(1-p_0)/n}}\sim N(0,1) \]

取$N(0,1)$的$1-\alpha$分位数为$u_{1-\alpha}$，则 $u\leq u_{1-\alpha}$时，认为总体的废品率$p$小于最低的废品率要求$p_0$

例 2. 甲方向乙方供货，约定废品率不超过3%，预先约定显著性水平$\alpha=0.05$。现在从一批产品中抽取100件，发现有5件废品，问乙方是否应该接受这批产品？

解. 知$\bar x=5/100$，$p_0=0.03$，$n=100$，可得

\[u=\frac{\bar x-p_0}{\sqrt{p_0(1-p_0)/n}}=1.1724 \]

而$u_{1-0.05}=1.65$。因此乙方应该接受这批产品。

Python来做假设检验

  • 正态性检验 
    scipy.stats.kstest(a_vector_like_data, 'norm') 
# 特点是比较严格，基于的原理是CDF，理论上可以检验任何分布。
scipy.stats.shapiro(a_vector_like_data)
# 专门用来检验正态分布。
scipy.stats.normaltest(a_vector_like_data)
# 原理是基于数据的skewness和kurtosis。
scipy.stats.anderson(a_vector_like_data, dist='norm')
# 是ks检验的正态检验加强版。
scipy.stats.jarque_bera(x)
# Perform the Jarque-Bera goodness of fit test on sample data.
  • 检验方差是否齐 
    对数据有正态性要求
scipy.stats.bartlett(a, b)
Levene检验,在数据非正态的情况下，精度比Bartlett检验好，可调中间值的度量
scipy.stats.levene(a, b, center = 'trimmed')
Fligner-Killeen检验,非参检验，不依赖于分布
scipy.stats.fligner(a, b, center='mean')
  • 两组数之间的比较
    1. 参数方法 
      # 独立两样本t检验
scipy.stats.ttest_ind(a, b, equal_var=True, nan_policy='omit')
# 成对两样本t检验
scipy.stats.ttest_rel(a, b, equal_var=True, nan_policy='omit')
# 通过基本统计量来做独立两样本检验
scipy.stats.ttest_ind_from_stats(20.06, 2.902, 50, 13.26, 1.977, 50, equal_var=False)
    2. 非参数方法 
      # wilcox秩序和检验，n < 20时独立样本效果比较好
scipy.stats.ranksums(a, b)
# Mann-Whitney U检验, n > 20时独立样本，比wilcox秩序和检验更稳健
scipy.stats.mannwhitneyu(a, b)
# Wilcox检验，成对数据
scipy.stats.wilcoxn(a, b, zero_method='wilcox', correction=False)
  • 多组数之间的比较 
    # 参数方法（1-way anova）
scipy.stats.f_oneway(a, b, c, ...)
# 非参数方法（Kruskal-Wallis H方法）
scipy.stats.kruskal(a, b, c,..., nan_policy='omit')
  • 相关性 相关性可以做简单的特征工程（特征筛选）来做监督学习以及作为相似度（1 - 距离）来做非监督学习。 
    # 参数（Pearson相关系数）
scipy.stats.pearsonr(a, b)
# 非参数（Spearman相关系数）
scipy.stats.spearmanr(a, b)
# 二元值和连续值之间的关系（Point-biserial相关系数）
scipy.stats.pointbiserialr(a, b)
#分参数的Kendall's Tau, 理论上是检验两个变量是否具有单调关系
scipy,stats.kendalltau(a, b, initial_lexsort=None, nan_policy='omit')

Python 

statsmodels.stats.stattools.jarque_bera

https://segmentfault.com/a/1190000007626742

Matlab

  • 单个总体 $N (\mu,\sigma) $ 均值 $\mu$ 的检验
    1. Z 检验法 
      [h,p,ci]=ztest(x,mu,sigma,alpha,tail)
    2. t 检验法 
      [h,p,ci]=ttest(x,mu,alpha,tail)
  • 两个正态总体均值差的检验( t 检验) 
    //还可以用 t 检验法检验具有相同方差的 2 个正态总体均值差的假设。
[h,p,ci]=ttest2(x,y,alpha,tail)

Matlab

  • 其它非参数检验 
    //Wilcoxon秩和检验
[p,h]=ranksum(x,y,alpha)
  • 中位数检验 
    //signrank函数,signrank Wilcoxon符号秩检验
[p,h]=signrank(x,y,alpha)
//signtest 符号检验
[p,h]= signtest(x,y,alpha)

例 3. 某车间用一台包装机包装糖果。包得的袋装糖重是一个随机变量,它服从正态分布。当机器正常时,其均值为 0.5 公斤,标准差为 0.015 公斤。某日开工后为检验包装机是否正常,随机地抽取它所包装的糖 9 袋,称得净重为(公斤) 

0.497, 0.506, 0.518, 0.524, 0.498, 0.511, 0.520, 0.515, 0.512

问机器是否正常?

例 4. 某种电子元件的寿命 $x$ (以小时计)服从正态分布, $\mu$, $\sigma$ 均未知.现得 16 只 元件的寿命如下: 

159, 280, 101, 212, 224, 379, 179, 264, 222, 362, 168, 250, 149, 260, 485, 170

问是否有理由认为元件的平均寿命大于 225(小时)?

例3 : matlab 

x=[0.497 0.506 0.518 0.524 0.498...
0.511 0.520 0.515 0.512];
[h,p,ci]=ztest(x,0.5,0.015) // 以mu=0.5, sigma=0.015来做假设，验证数据

求得 h=1,p=0.0248。 $p<0.05$，说明可拒绝原假设，即认为这天包装机工作不正常。

例4 : matlab 

x=[159 280 101 212 224 379 179 264 ...
222 362 168 250 149 260 485 170];
[h,p,ci]=ttest(x,225,0.05,1) // 以 mu=225 来做假设，验证数据

求得 h=0,p=0.2570。 $p>0.05$，说明不能 拒绝原假设，即认为元件的平均寿命不大于 225 小时。

Python 

import numpy as np
from scipy import stats as st

data4=np.array([159, 280, 101, 212, 224, 379, 179, 264, 222, 362, 168, 250, 149, 260, 485, 170.0])

s,p=st.kstest(data4,'norm')
print('kstest:',s,p)
# 使用scipy包中的ttest_1samp,计算单独样本t检验。该函数返回的第一个值为t值，第二个值为双尾检验的p值。
t,p=st.ttest_1samp(data4,225.0)
 # 建立判断标准检验证据是否有效，给定的判断标准即显著水平
 # 当p<=5%*2时，拒绝零假设，接受备选假设，即符合要求
 #   p>5%*2时，接受零假设
if(p<=0.01): print("样本满足正态分布假定")
else: print("样本*不*满足正态分布假定")
ttest: 0.668517696746 0.513960143175

Python 的 scipy.stats 库没有ztest函数，可以自己写 

import numpy as np
from scipy import stats as st

data=np.array([0.497, 0.506, 0.518, 0.524, 0.498, 0.511, 0.520, 0.515, 0.512])

# z-test for data with N(0.5, 0.015)
def ztest(data, mu, sigma):
    barx=data.mean()
    # print(np.sqrt(len(data)))
    z=(barx-mu)/(sigma/np.sqrt(len(data)))
    # 将Z分数(Z值,标准分数)转换为Python中正态分布的p值
    # https://codeday.me/bug/20180228/136914.html
    p_values=st.norm.sf(abs(z))*2 # two-sided,
    #p2_val = st.norm.sf(st.norm.ppf(1-0.025))
    return {"z":z, "p":p_values, "mean":barx}
print("======= Z-test of data =========")
res = ztest(data, 0.5, 0.015)
print("Z=",res['z'], ", p=",res['p'])

得到结果 

Z= 2.24444444444 , p= 0.0248038196323
statsmodels.stats.weightstats.ztest(x1, x2=None, value=0, alternative='two-sided',
                                    usevar='pooled', ddof=1.0)

方差分析

例 5. (1993美国数模竞赛A) 有堆肥数据如下

mm7-ex-umm1993a

哪些因素与堆肥形成的速度有关？

实际中，许多问题涉及多因素的影响，一般要先在众多的因素中选出影响大的，以进一步做更细致的试验。用来判断一个因素的影响“是否大”的主要方法就是方差分析(Analysis of Variance，简称ANOVA)。

方差分析(Analysis of Variance，简称ANOVA)，又称“变异数分析”，是R.A.Fisher发明的，用于两个及两个以上样本均数差别的显著性检验。 由于各种因素的影响，研究所得的数据呈现波动状。造成波动的原因可分成两类，一是不可控的随机因素，另一是研究中施加的对结果形成影响的可控因素。方差分析是从观测变量的方差入手，研究诸多控制变量中哪些变量是对观测变量有显著影响的变量。

方差分析分为单因素方差分析和多因素分析。主要是单因素方差分析

单独考察某一个因素$A$对指标$X$的影响，可以采用单因素方差分析。

设因素$A$有$k$个水平，每一个水平下做了$n$次实验，得到数据

水平$h_i$： $X_{i1}$, $X_{i2}$, $\cdots$, $X_{in}$，样本均值$\displaystyle\bar X_i=\frac1n\sum_{j=1}^nX_{ij}$

在每一个水平下，认为$X_{i1}$,$\cdots$, $X_{in}$之间的差异与因素$A$无关，只是由随机误差所引起。通常，假定$X_{i1}$,$\cdots$, $X_{in}$服务正态分布$N(\mu_i,\sigma_i^2)$。通常$\mu_1$, $\cdots$, $\mu_k$不同，但是，如果因素$A$对指标$X$的影响不大，则可以认为$\mu_1=\cdots=\mu_k$。因此，给出假设检验：

\[H_0: \mu_1=\mu_2=\cdots=\mu_k \]

若返回$p$值小于0.05，则假设不成立，即因素$A$对指标$X$影响很大；否则，表示影响不大。

Python 

scipy.stats.f_oneway(x1,x2,x3,...)
from statsmodels.stats.anova import anova_lm

返回F值和p值。若显示p值小于0.05，则认为它们有显著性差异。

多因素方差分析

Python 

from statsmodels.stats.anova import anova_lm
from statsmodels.formula.api import ols
anova_result=anova_lm(ols(...),typ=)

p值小于0.05的，表示该因素有影响 

ref： https://www.cnblogs.com/tangxianwei/p/8359472.html

Matlab 统计工具箱中单因素方差分析的命令是 

anoval(x)

用 anova2 作双因素方差分析。命令为 

p=anova2(x,reps)

可以考虑的因素有很多：混合物的总量，混合物中各成份的比例，环境温度，微生物的数量，等等。

  1. 温度影响（考虑堆肥的时间）

    把时间分成7月分，8月份，4-5月份，共3个水平，有数据： 
    水平1： 28,27,27,26
水平2： 33,36,35,47
水平3： 49,49,49,49
    分析后有 
    F_onewayResult(statistic=36.08074534161494, pvalue=5.035194192575532e-05)
    可以看到p值小于0.05，因此，温度因素（堆肥的时间）对堆肥形成时间影响显著。 

from scipy import stats
import numpy as np
# 时间的影响
c_d1=np.array([28.0,27,27,26])
c_d2=np.array([33.0,36,35,47])
c_d3=np.array([49.0,49,49,49])
c_d4=np.array([28.0,27,27,26])
# 单因素方差分析
res=stats.f_oneway(c_d1,c_d2,c_d3)
print(res)
# 结果 pvalue=0.00005 < 0.05 ，即因素有显著性差异，表明因素有影响

  1. C/N比对堆肥形成时间的影响

    将C/N比分为15~16，23~24，25以上共3个水平，对应形成天数为 
    水平1：28,27,27,26,33,36,35,47
水平2： 49,49
水平3： 49,49
    分析后，有 
    F_onewayResult(statistic=9.425754884547072, pvalue=0.006198290861822317)
    p值小于0.05，可以认为有影响
# C/N比
print("----------C/N")
c_d1=np.array([28.0,27,27,26,33,36,35,47]) # 取值 15
c_d2=np.array([49.0,49])  # 取值 23
c_d3=np.array([49.0,49])  # 取值 25-26
# 单因素方差分析
res=stats.f_oneway(c_d1,c_d2,c_d3)
print(res)
# pvalue=0.0061 < 0.05 ，影响显著

  1. 混合物总重量对堆肥形成时间的影响

    将总重量分为5个水平：91～115，116～140，141～165，191～215，266～290 
    水平1： 27,33,49 ；水平2： 28,35,49,49,49
水平3： 36,47； 水平4： 27； 水平5： 26
    可以得到结果 
    F_onewayResult(statistic=0.9917879914674375, pvalue=0.4786196984781391)
    显示p>0.05，因此没有影响
c_d1=np.array([27.0,33,49]) # 91~115
c_d2=np.array([28.0,49,49,49]) # 116~140
c_d3=np.array([36.0,47]) # 141~165
c_d4=np.array([27.0]) # 191~215
c_d5=np.array([26.0]) # 266~290
# 单因素方差分析
res=stats.f_oneway(c_d1,c_d2,c_d3,c_d4,c_d5)
print(res)
# pvalue=0.47 > 0.05 ，没有影响

例 6. 为考察 5 名工人的劳动生产率是否相同,记录了每人 4 天的产量。你能从这些数据推断出他们的生产率有无显著差别吗?

工人 $A_1$ $A_2$ $A_3$ $A_4$ $A_5$
第1天 256 254 250 248 236
第2天 242 330 277 280 252
第3天 280 290 230 305 220
第4天 298 295 302 289 252 
c_d1=np.array([256.0, 242, 280, 298]) # dose 1
c_d2=np.array([254.0, 330, 290, 295]) # dose 2
c_d3=np.array([250.0, 277, 230, 302]) # dose 3
c_d4=np.array([248.0, 280, 305, 289]) # dose 2
c_d5=np.array([236.0, 252, 220, 252]) # dose 3
## 单因素方差分析
res=stats.f_oneway(c_d1,c_d2,c_d3,c_d4,c_d5)
print(res)

可以得到 

F_onewayResult(statistic=2.2617412494461675, pvalue=0.1109133645884241)

p=0.11>0.05，所以没有差异。

回归分析

例 7. 某金融公司打算新开一类金融产品，现有9个金融产品的数据，包括用户购买金融产品的综合年化利率，以及公司收取用户的佣金（手续费）；如下表所示，产品利率为11％，佣金为50，我们需要预测这款金融产品的销售额

产品 1 2 3 4 56789 10
年化利率 9 7 756786911
用户佣金 7530203002150206050
销售额 500370375270360379440300510？

解. 销售额与年化利率和用户佣金之间的关系是什么？

回归分析（regression analysis)是研究随机变量之间的关系。

具体地说,回归分析在一组数据的基础上研究这样几个问题:

  1. 建立因变量 $y$ 与自变量 $x_1$, $\cdots$, $x_m$ 之间的回归模型(经验公式);
  2. 对回归模型的可信度进行检验;
  3. 判断每个自变量 $x_i$ 对 $y$ 的影响是否显著;
  4. 诊断回归模型是否适合这组数据;
  5. 利用回归模型对 $y$ 进行预报或控制。

如果模型中只含有一个回归变量，称为一元回归模型，否则，称为多元回归模型。

  • 曲线拟合问题的特点是,根据得到的若干有关变量的一组数据,寻找因变量与(一个或几个)自变量之间的一个函数,使这个函数对那组数据拟合得最好。
  • 通常,函数的形式可以由经验、先验知识或对数据的直观观察决定,要作的工作是由数据用最小二乘法计算函数中的待定系数。 从计算的角度看,问题似乎已经完全解决了,还有进一步研究的必要吗?

看一个虚拟的例子：

mm7-linear-fit

可以看到，$z$与$x$的模型是没有什么价值的。

  • 用统计方法对结果进行定理研究是回归分析的任务。简单地说,回归分析就是对拟合问题作的统计分析。
  • 从数理统计的观点看,这里涉及的都是随机变量,我们根据一个样本计算出的那些系数,只是它们的一个(点)估计,应该对它们作区间估计或假设检验,如果置信区间太大,甚至包含了零点,那么系数的估计值是没有多大意义的。
  • 另外也可以用方差分析方法对模型的误差进行分析,对拟合的优劣给出评价。

线性回归

一元线性回归模型可以表示为

\[y=\beta_0+\beta_1x+\epsilon \]

其中$\beta_0$, $\beta_1$待定，称为回归系数。其基本假设是，对于不同的$x$随机误差$\epsilon$是相互独立的随机变量，且

\[\epsilon \sim N(0,\sigma^2) \]

多元线性回归的模型为

\[\begin{cases} y=\beta_0+\beta_1x_1+\cdots+\beta_mx_m+\epsilon \\ \epsilon \sim N(0,\sigma^2) \end{cases} \]

一般地,回归方程的假设检验包括两个方面:

  • 一个检验是回归参数的检验-$t$检验
    1. 这个检验的是每一个自变量对因变量的影响程度是否显著。
    2. 若系数$\beta_i$, $i>0$的这个置信区间包含了0，则意味着$\beta_i=0$是可以接受的，就是说$y$依赖于$x_i$的那一项毫无意义。
    3. 对于常系数$\beta_0$的区间估计要复杂些，意义也不大。
  • 一个是对模型的检验-$F$检验
    1. 检验自变量与因变量之间的关系能否用一个线性模型来表示
    2. 可以得到决定系数$R^2$参数，用来衡量模型的可靠性
  • 从逻辑上说,一般常在 F 检验通过后,再进一步进行 t 检验。

由最小二乘法，可以得到

\[\hat y=\hat\beta_0+\hat\beta_1x_1+\cdots+\hat\beta_mx_m \]

则对样本数据，有

\[Q=\sum_{i=1}^n e_i^2=\sum_{i=1}^n(y_i-\hat y_i)^2 \]

称为残差平方和(或剩余平方和)。

总平方和$SST=\sum_{i=1}^n(y_i-\bar y)^2$可以分解为

\[SST=Q+U, U=\sum_{i=1}^n(\hat y_i-\bar y)^2 \]

$U$称为回归平方和，反映自变量对 $y$ 的影响。

给出假设

\[H_0: \beta_j=0, j=1,2,\cdots,m \]

（也就是说，$y$与变量$x_i$的线性关系不明显。）在$H_0$假设下，有

\[F=\frac{U/m}{Q/(n-m-1)}\sim F(m,n-m-1) \]

在显著性水平 $\alpha$ 下，若$F<F_{\alpha}(m,n-m-1)$，接受$H_0$；否则，拒绝$H_0$。

另外，用

\[R^2=\frac{U}{SST} \]

来衡量 $y$ 与 $x_i$的相关程度。 $R$越大，关系越密切。通常要$R>0.8$甚至$0.9$。

在上面$H_0$被拒绝时，也就是$\beta_j$不全为零。但不排除其中若干个等于零。进一步作如下 $m + 1$ 个检验

\[H^{(j)}_0: \beta_j=0, j=0,1,\cdots,m \]

系数$\hat\beta_j$满足$\hat\beta_j\sim N(\beta_j,\sigma^2 c_{jj})$，$c_{jj}$是$(X^TX)^{-1}$中的第$(j,j)$元素，$X^TX$是最小二乘时产生的矩阵。用剩余方差$s^2$代替$\sigma^2$

\[s^2=\frac{Q}{n-m-1} \]

可以得到， 当$H_0^{(j)}$成立时，

\[t_j=\frac{\hat\beta_j/\sqrt{c_{jj}}}{\sqrt{Q/(n-m-1)}}\sim t(n-m-1) \]

对于给定的$\alpha$，若$|t_j|<t_{\frac{\alpha}2}(n-m-1)$，接受$H^{(j)}_0$；否则，拒绝。

同时，式子也可以给出$\hat\beta_j$的置信区间

\[[\hat\beta_j-t_{\frac{\alpha}2}(n-m-1), \hat\beta_j+t_{\frac{\alpha}2}(n-m-1)] \]

Matlab 统计工具箱用命令 regress 实现多元线性回归,用的方法是最小二乘法,用 法是: 

b=regress(Y,X)

其中 Y,X 为数据,b 为回归系数估计值 。 

[b,bint,r,rint,stats]=regress(Y,X,alpha)

这里 Y,X 同上,alpha 为显著性水平(缺省时设定为 0.05),b,bint 为回归系数估计值和 它们的置信区间,r,rint 为残差(向量)及其置信区间,stats 是用于检验回归模型的统计量

Python 

# 简单的线性回归分析，不能用它拟合一般的线性模型，或者是用它来进行多变量回归分析
scipy.stats.linregress(x, y=None)
#
optimize.curve_fit( )
# statsmodels
import statsmodels.api as sm
model = sm.OLS(y, x).fit()
# sklearn
from sklearn.linear_model import LinearRegression
lrModel = LinearRegression()
lrModel.fit(x,y)

Mathemtica 

model = LinearModelFit[data, x, x]
model["BestFit"] (* 提取模型的函数形式 *)
Show[ListPlot[data], Plot[model["BestFit"], {x, 0, 30}]] (* 绘制模型的函数形式 *)
model["ParameterTable"] (* 获取有关参数估计的信息：P-value等 *)
model["ParameterConfidenceIntervalTable"]  (* 获取参数的区间估计 *)

例 8. 某金融公司打算新开一类金融产品，现有9个金融产品的数据，包括用户购买金融产品的综合年化利率，以及公司收取用户的佣金（手续费）；如下表所示，产品利率为11％，佣金为50，我们需要预测这款金融产品的销售额

产品 1 2 3 4 56789 10
年化利率 9 7 756786911
用户佣金 7530203002150206050
销售额 500370375270360379440300510？

用sklearn.linear_model的结果 

决定系数（贡献率）： 0.960315344209031
预测销售额是：  [[627.33389606]]
回归系数是：beta0=[-44.5394503], beta_i=[[62.35235754 -0.28005173]]

可以看到决定系数离1很近，表明模型的有效性很大。

用statsmodels.api.OLS的结果 

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.960
Model:                            OLS   Adj. R-squared:                  0.947
Method:                 Least Squares   F-statistic:                     72.60
Date:                Thu, 11 Apr 2019   Prob (F-statistic):           6.25e-05
Time:                        10:53:36   Log-Likelihood:                -37.326
No. Observations:                   9   AIC:                             80.65
Df Residuals:                       6   BIC:                             81.24
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        -44.5395     46.869     -0.950      0.379    -159.224      70.145
x1            62.3524      8.253      7.555      0.000      42.158      82.547
x2            -0.2801      0.482     -0.581      0.583      -1.460       0.900
==============================================================================
Omnibus:                        0.629   Durbin-Watson:                   2.077
Prob(Omnibus):                  0.730   Jarque-Bera (JB):                0.302
Skew:                           0.392   Prob(JB):                        0.860
Kurtosis:                       2.562   Cond. No.                         312.
==============================================================================
Prediction by statsmodels: 
[627.33389606]
import numpy as np

x = np.array([[9.0,75.0],[7,30],[7,20],[5,30],[6,0],[7,21],[8,50],[6,20],[9,60]]) # data[["百分比利率","抽取用户佣金"]]
y = np.array([[500.0,], [370,], [375,], [270,], [360,], [379,], [440,], [300,],[510,]]) #data[["金融产品销售额"]]

#建模
from sklearn.linear_model import LinearRegression
lrModel = LinearRegression()
#训练模型
lrModel.fit(x,y)
## 查看R^2，它在(0,1)之间，越大说明回归模型越有效
print("决定系数（贡献率）：")
print(lrModel.score(x,y))
#预测
amount=lrModel.predict(np.array([[11,50]]))
print("预测销售额是： ", amount)
print("回归系数是：")
#查看截距  \beta_0
print(lrModel.intercept_)
#查看参数 \beta_1,\beta_2,\cdots
print(lrModel.coef_)
import statsmodels.api as sm
#print(x)
x = sm.add_constant(x) ## let's add an intercept (beta_0) to our model
#print(x)
model = sm.OLS(y, x).fit()
# Print out the statistics
print(model.summary())
print("Prediction by statsmodels: ")
predictions = model.predict(np.array([[1.0,11,50]])) # make the predictions by the model
print(predictions)

一般的多元线性回归模型

\[\begin{cases} r(y)=\beta_0+\beta_1r_1(x)+\cdots+\beta_mr_m(x)+\epsilon \\ \epsilon \sim N(0,\sigma^2) \end{cases} \]

如

\[y=\beta_0+\beta_1x^2, y=\beta_0+\beta_1e^{x_1}+\beta_2\frac1{x_2} \]

对$y=\beta_0e^{\beta_1x}$，可以做

\[\ln y=\beta_0+\beta_1x \]

的线性回归

多元二项式回归

  1. linear(线性): $y=\beta_0+\beta_1x_1+\cdots+\beta_mx_m$
  2. purequadratic(纯二次):
    \[\displaystyle y =\beta_0+\beta_1x_1+\cdots+\beta_mx_m+\sum_{j=1}^m\beta_{jj}x_j^2 \]
  3. interaction(交叉):
    \[\displaystyle y =\beta_0+\beta_1x_1+\cdots+\beta_mx_m+\sum_{1\leq j<k\leq m}\beta_{jk}x_jx_k \]
  4. quadratic(完全二次):
    \[\displaystyle y =\beta_0+\beta_1x_1+\cdots+\beta_mx_m+\sum_{1\leq j\leq k\leq m}\beta_{jk}x_jx_k \]

Python 

# 单变量高次多项式的拟合
numpy.polyfit(x,y,deg) # y=a_0+a_1 x+...+a_deg x^deg

# 多变量高次多项式的线性回归
from sklearn.linear_model import LinearRegression#导入线性回归模型
from sklearn.preprocessing import PolynomialFeatures# 导入多项式回归模型
quadratic_featurizer = PolynomialFeatures(degree = 2)#实例化一个二次多项式
x_train_quadratic = quadratic_featurizer.fit_transform(x)#用二次多项式多样本x做变换
# 把训练好X值的多项式特征实例应用到一系列点上,形成矩阵
xx_quadratic = quadratic_featurizer.transform(xx.reshape(xx.shape[0], 1))

regressor_quadratic = LinearRegression()
regressor_quadratic.fit(x_train_quadratic, y) # 线性回归

非线性回归

非线性回归是指因变量 $y$ 对回归系数 $\beta_1$,$\cdots$, $\beta_m$ (而不是自变量)是非线性的。 如

\[y=\frac{\beta_4x_2-\frac{x_3}\beta_5}{\beta_1x_1+\beta_2x_2+\beta_3x_3} \]
  • 可以用Gauss-Newton法来求相应的最小二乘问题
  • Matlab统计工具箱中的nlinfit计算回归系数,用nlparci计算回归系数的置信区间,用nlpredci计算预测值及其置信区间,
  • Python中，可以用scipy.optimize.curve_fit(...)来求回归系数。

例 9. (软件开发人员的薪金) 研究薪金与人员资历、管理责任、教育程度的关系

编号 薪金 资历 管理 中学教育 大学教育
1 13876 1 1 1 0
2 11608 1 0 0 0
3 18701 1 1 0 0
… …
46 19346 20 0 1 0

解. 薪金为$y$，资历为$x_1$，管理为$x_2$，教育程度由中学教育与大学教育共同定义：中学用$(1,0)$，大学用$(0,1)$，研究生用$(0,0)$。先建立线性模型

\[y=a_0+a_1x_1+a_2x_2+a_3x_3+a_4x_4 \]

用Python的statsmodels的计算结果 

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.957
Model:                            OLS   Adj. R-squared:                  0.952
Method:                 Least Squares   F-statistic:                     226.4
Date:                Sun, 14 Apr 2019   Prob (F-statistic):           2.31e-27
Time:                        11:49:24   Log-Likelihood:                -381.66
No. Observations:                  46   AIC:                             773.3
Df Residuals:                      41   BIC:                             782.5
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const       1.103e+04    383.492     28.769      0.000    1.03e+04    1.18e+04
x1           546.1276     30.541     17.882      0.000     484.449     607.807
x2          6882.5329    314.145     21.909      0.000    6248.105    7516.961
x3         -2994.1783    412.049     -7.267      0.000   -3826.327   -2162.029
x4           147.7380    387.938      0.381      0.705    -635.718     931.194
==============================================================================
import numpy as np
a=np.loadtxt("data/data0902.txt")
# 得到数据（第2列）
y=a[:,1]
x=a[:,[2,3,5,6]]
# ============ 用statsmodels做线性回归 -------------
import statsmodels.api as sm
#print(x)
x = sm.add_constant(x) ## let's add an intercept (beta_0) to our model
#print(x)
model = sm.OLS(y, x).fit()
# Print out the statistics
print(model.summary())

结果分析

  1. R-squared: 0.957表明薪水95.7%可以由模型来计算，因而整体来说，模型是可用的
  2. x1: 546.1276表明每涨一年资历，工资涨546
  3. x2: 6882.5329表明，管理人员的薪水比非管理人员高6883
  4. x3: -2994表明，中学程度的薪水比研究生少 2994
  5. x4: 148表明，大学程度的薪水比研究生多148。但是这个系数的置信区间是-635.718, 931.194，包含了0，所以这个系数不太可靠

进一步分析: x4的系数的置信区间包含了0，模型还存在缺点。用残差分析法。（残差是薪水的实际值$y$与模型的估计值$\hat y$之差）。

把影响因素分为资历，管理-教育2类。管理-教育按管理+教育（中学为1，大学为2，研究生为3）给出值。

先看资历：残差图。大概分为3个水平，这是因为混合了管理-教育。

mm7-ex-salary-err1

# 计算残差 \epsilon
err=y-model.fittedvalues

import matplotlib.pyplot as plt
# 画出 err-资历 图
plt.plot(x[:,1],err,'+')
plt.show()
# 画出 err- 管理-学历组合图
plt.xlim(0.5,6.7)
plt.plot(a[:,7],err,"+")
plt.show()

再看教育-管理：残差图。残差全为正或全为负，表明管理-教育组合在模型中处理不当

mm7-ex-salary-err2

更好的模型： 增加管理$x_2$与教育程序$x_3$, $x_4$的交互项，

\[y=a_0+a_1x_1+a_2x_2+a_3x_3+a_4x_4+a_5x_2x_3+a_6x_2x_4 \]

可以得到结果 R-squared: 0.999和 

                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        1.12e+04     78.853    142.085      0.000     1.1e+04    1.14e+04
x1           496.8639      5.551     89.502      0.000     485.635     508.093
x2          7047.9997    102.313     68.887      0.000    6841.052    7254.948
x3         -1726.5042    105.050    -16.435      0.000   -1938.989   -1514.020
x4          -348.3925     97.305     -3.580      0.001    -545.211    -151.574
x5         -3070.5962    148.929    -20.618      0.000   -3371.833   -2769.360
x6          1835.9676    130.814     14.035      0.000    1571.370    2100.565
==============================================================================

可以看到，$R^2$值和系数的置信区间都比前一个模型好

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.999
Model:                            OLS   Adj. R-squared:                  0.999
Method:                 Least Squares   F-statistic:                     5545.
Date:                Sun, 14 Apr 2019   Prob (F-statistic):           1.51e-55
Time:                        15:43:29   Log-Likelihood:                -298.62
No. Observations:                  46   AIC:                             611.2
Df Residuals:                      39   BIC:                             624.0
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        1.12e+04     78.853    142.085      0.000     1.1e+04    1.14e+04
x1           496.8639      5.551     89.502      0.000     485.635     508.093
x2          7047.9997    102.313     68.887      0.000    6841.052    7254.948
x3         -1726.5042    105.050    -16.435      0.000   -1938.989   -1514.020
x4          -348.3925     97.305     -3.580      0.001    -545.211    -151.574
x5         -3070.5962    148.929    -20.618      0.000   -3371.833   -2769.360
x6          1835.9676    130.814     14.035      0.000    1571.370    2100.565
==============================================================================
# --------- 增加交互项 x2*x3, x2*x4
#print(x)
x0=x[:,2]*x[:,3]
x1=x[:,2]*x[:,4]
x=np.column_stack((x,x0.T))
x=np.column_stack((x,x1.T))

model2 = sm.OLS(y, x).fit()
# Print out the statistics
print(model2.summary())

资历：残差图。

mm7-ex-salary-m2-err1

教育-管理：残差图。

mm7-ex-salary-m2-err2

两幅图中，都可以看到一个异常点：10年资历的某个人（查数据，为33号数据）。他的实际薪水明显低于有类似经历的其他人的薪水。应该把这个异常数据去掉。 

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       1.000
Model:                            OLS   Adj. R-squared:                  1.000
Method:                 Least Squares   F-statistic:                 3.670e+04
Date:                Sun, 14 Apr 2019   Prob (F-statistic):           6.65e-70
Time:                        16:09:26   Log-Likelihood:                -248.54
No. Observations:                  45   AIC:                             511.1
Df Residuals:                      38   BIC:                             523.7
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        1.12e+04     29.995    373.401      0.000    1.11e+04    1.13e+04
x1           498.2944      2.114    235.743      0.000     494.015     502.573
x2          7041.1690     38.920    180.914      0.000    6962.380    7119.958
x3         -1737.0900     39.965    -43.466      0.000   -1817.994   -1656.186
x4          -356.3558     37.016     -9.627      0.000    -431.291    -281.420
x5         -3056.3268     56.657    -53.945      0.000   -3171.022   -2941.632
x6          1996.9828     50.871     39.256      0.000    1893.999    2099.966

资历：残差图。

mm7-ex-salary-m3-err1

教育-管理：残差图。

mm7-ex-salary-m3-err2

  1. 高中学历的起薪是多少？ 9463
  2. 高中学历的管理岗位的起薪是多少？ 13448
  3. 大学学历的起薪是多少？ 10844
  4. 大学学历的管理岗位的起薪是多少？ 19882
  5. 研究生学历的起薪是多少？ 11200
  6. 研究生学历的管理岗位的起薪是多少？ 18241

Logistic回归模型

例 10. 为了更好地说明冠心病发病率与年龄的关系，医学家们对100名不同年龄的人进行观察，得到这100名被观察者的年龄及他们是否患冠心病的数据。下面分析冠心病发病率与年龄的关系，并进行预测

序号 1 2 … 100
年龄 20 23 … 69
冠心病 0 0 … 1

数据的散点图

mm7-ex-heart-ill-1

直接回归是不可能的。先做预处理

数据预处理： 按年龄段进行分组

年龄段 年龄中点 人数 患病人数 患病比例
20~29 24.5 10 1 0.1
30~34 32 15 2 0.13
35~39 37 12 3 0.25
40~44 42 15 5 0.33
45~49 47 13 6 0.46
50~54 52 8 5 0.63
55~59 57 17 13 0.76
60~69 64.5 10 8 0.8

可以做出1次和3次的回归曲线

\[y=\beta_0+\beta_1 x+\beta_2 x^2+\beta_3 x^3+\epsilon \]

mm7-ex-heart-lin mm7-ex-heart-poly

问题：

  1. $y$值可能超过[0,1]
  2. $\epsilon$不具有正态性
import statsmodels.api as sm

X=sm.add_constant(Age) ## let's add an intercept (beta_0) to our model
y=proportion

# 对年龄段中点-患病比例做 linear regression
lin_mod = sm.OLS(y, X)
lin_res = lin_mod.fit()
# polynomial regression
degree=3
weights = np.polyfit(Age.ravel(), y, degree)
print(weights)
poly_mod = np.poly1d(weights)

当$Y$是一个二分类（或多分类）变量，而不是连续变量时，普通的回归分析是不适合的，需要使用logistic回归模型。

用$\pi(x)$表示年龄为$x$的患病概率，

\[\pi(x)=P(Y=1|x) \]

为了得到$\pi(x)$与$x$函数关系，同时注意到$\pi(x)$在$[0,1]$间取值，使用Logistic模型

\[\pi(x)=\frac1{1-e^{\beta_0+\beta_1x}} \]

即

\[\ln\frac{\pi(x)}{1-\pi(x)}=\beta_0+\beta_1x \]

属于广义线性模型。（Generalized Linear Model）

记

\[logit(\pi(x))=\ln\frac{\pi(x)}{1-\pi(x)} \]

称为logit模型或logistic回归模型（二项式逻辑回归模型）（分类评定模型）（logistic regression）

Matlab 

b=glmfit(x,y,'binomial','logit')

Mathematica 

logit = LogitModelFit[data, x, x]
Normal[logit] (* 查看模型的函数形式 *)
Show[ListPlot[data], Plot[logit[x], {x, 1, 3}], PlotRange -> All, AxesOrigin -> {.9, -.1}] (*绘制数据点和模型*)
(* GeneralizedLinearModelFit 中一个缺省的 "Binomial" 模型等价于 LogitModelFit 的模型 *)
GeneralizedLinearModelFit[data, x, x, ExponentialFamily -> "Binomial"] // Normal

Python 

import statsmodels.api as sm
logit_mod = sm.Logit(y,X)
logit_res=logit_mod.fit()
#logit_res.predict([1, t])
#logit_res.summary()
from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
X, y = load_iris(return_X_y=True)
clf = LogisticRegression(random_state=0, solver='lbfgs',
                         multi_class='multinomial').fit(X, y)
clf.predict(X[:2, :])
# array([0, 0])
clf.predict_proba(X[:2, :]) 
array([[9.8...e-01, 1.8...e-02, 1.4...e-08],
       [9.7...e-01, 2.8...e-02, ...e-08]])
clf.score(X, y)
# 0.97...
# https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression

结果

mm7-ex-heart-logit mm7-ex-heart-3

                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -4.8981      3.749     -1.306      0.191     -12.246       2.450
x1             0.1019      0.079      1.287      0.198      -0.053       0.257
==============================================================================

记$odds(x)$为年龄$x$的人患病与不患病的概率之比，则

\[odds(x)=\frac{\pi(x)}{1-\pi(x)}=e^{\beta_0+\beta_1x} \]

各年龄的发病概率与发生比

\[\begin{aligned} \pi(20)= [0.05416873] [0.05727103] \\ \pi(30)= [0.13694377] [0.15867305] \\ \pi(40)= [0.30536927] [0.43961382] \\ \pi(50)= [0.54913894] [1.21797818] \\ \pi(60)= [0.77140173] [3.37448636] \\ \end{aligned} \]

$\pi(48)= [0.49834553]$，$\pi(49)= [0.52380378]$。 在48岁后，患病的概率会大于不患病的概率。

除了使用logit模型来处理这种0-1变量模型外，还可以用另一种广义线性模型probit模型，其形式为

\[probit(\pi(x))=\Phi^{-1}(\pi(x))=\beta_0+\beta_1x \]

其中$\Phi$是正态概率分布函数，它也是$S$形曲线。

Matlab 

[probitCoef,devp,statsp] = glmfit(Age,[Chd Total],'binomial','probit');  % probit 模型求解

Mathematica 

probit = ProbitModelFit[data, x, x]
(* ProbitModelFit 等价于有 "ProbitLink" 的一个 "Binomial" 模型 *)
GeneralizedLinearModelFit[data, x, x, ExponentialFamily -> "Binomial",
LinkFunction -> "ProbitLink"] // Normal

Python 

import statsmodels.api as sm

probit_mod = sm.Probit(y,X)
probit_res = probit_mod.fit()
probit_margeff = probit_res.get_margeff()
print('Parameters: ', probit_res.params)
print('Marginal effects: ')
print(probit_margeff.summary())

目录

Stopping stepwise: Why stepwise selection is bad and what you should use instead

https://towardsdatascience.com/stopping-stepwise-why-stepwise-selection-is-bad-and-what-you-should-use-instead-90818b3f52df

LASSO and LAR offer much better alternatives than stepwise as a starting point for further analysis.

Python offers LASSO through sklearn: 

>>> from sklearn import linear_model
>>> clf = linear_model.Lasso(alpha=0.1)
>>> clf.fit([[0,0], [1, 1], [2, 2]], [0, 1, 2])
Lasso(alpha=0.1, copy_X=True, fit_intercept=True, max_iter=1000,
   normalize=False, positive=False, precompute=False, random_state=None,
   selection='cyclic', tol=0.0001, warm_start=False)
>>> print(clf.coef_)
[ 0.85  0.  ]
>>> print(clf.intercept_)
0.15
GLS.fit_regularized(method='elastic_net', alpha=0.0, L1_wt=1.0, start_params=None, profile_scale=False, refit=False, **kwargs)
# Return a regularized fit to a linear regression model.
# method (string) – ‘elastic_net’ and ‘sqrt_lasso’ are currently implemented.
# L1_wt (scalar) – The fraction of the penalty given to the L1 penalty term. Must be between 0 and 1 (inclusive). If 0, the fit is a ridge fit, if 1 it is a lasso fit.