进行参数估计和假设检验时，通常总是假定总体服从正态分布，虽然在许多情况下这个假定是合理的，但是当要以此为前提进行重要的参数估计或假设检验，或者人们对它有较大怀疑的时候，就确有必要对这个假设进行检验，进行总体正态性检验的方法有很多种，以下针对MATLAB统计工具箱中提供的程序，简单介绍几种方法。
1)Jarque-Bera检验
利用正态分布的偏度g1和峰度g2，构造一个包含g1，g2的分布统计量(自由度n=2)，对于显著性水平，当分布统计量小于分布的分位数时，接受H0：总体服从正态分布；否则拒绝H0，即总体不服从正态分布。这个检验适用于大样本，当样本容量n较小时需慎用。Matlab命令：h
=jbtest(x)，[h,p,jbstat,cv] =jbtest(x,alpha)。
2)Kolmogorov-Smirnov检验
通过样本的经验分布函数与给定分布函数的比较，推断该样本是否来自给定分布函数的总体。容量n的样本的经验分布函数记为Fn(x)，可由样本中小于x的数据所占的比例得到，给定分布函数记为G(x)，构造的统计量为，即两个分布函数之差的最大值，对于假设H0：总体服从给定的分布G(x)，及给定的，根据Dn的极限分布(n??时的分布)确定统计量关于是否接受H0的数量界限。
因为这个检验需要给定G(x)，所以当用于正态性检验时只能做标准正态检验，即H0：总体服从标准正态分布。Matlab命令：h
=kstest(x)。
3)Lilliefors检验
它将Kolmogorov-Smirnov检验改进用于一般的正态性检验，即H0：总体服从正态分布，其中由样本均值和方差估计。Matlab命令：
h =lillietest(x)，[h,p,lstat,cv]=lillietest(x,alpha)。
4)另外还有一种方法：首先对于数据进行标准化：Z =
ZSCORE(X)，然后在进行2)的Kolmogorov-Smirnov检验，检验是否为标准正态分布，类似于对于方法2)的改进。

文章目录
九、独立性检验和正态性检验
1.独立性检验
2.一元数据正态性检验
3.多元数据的正态性检验
回顾总结

九、独立性检验和正态性检验
1.独立性检验
独立性检验，指的是将一个多元总体$X\sim N_p(\mu,\Sigma)$划分成$k$个部分，探究每个部分之间是否独立的问题，这样做的好处是显而易见的，如果一个总体$X$可以划分成多个独立的部分，那么只需要对每一个部分分开讨论即可，无疑降低了运算量。在多元统计中，可以视为有如下分解：

$X=\begin{bmatrix}
X^{(1)}  \\ \vdots \\ X^{(k)}
\end{bmatrix},
\mu=\begin{bmatrix}
\mu^{(1)}  \\ \vdots \\ \mu^{(k)}
\end{bmatrix}, 
\Sigma=\begin{bmatrix}
\Sigma_{11} & \cdots & \Sigma_{1k} \\
\vdots & & \vdots \\
\Sigma_{k1} & \cdots & \Sigma_{kk}
\end{bmatrix}.$

每一个分向量$X^{(t)}$都是$p_t$维的，对应的$\mu^{(t)}$也是$p_t$维的，$\Sigma_{tt}$是$p_t\times p_t$的。在多元正态分布的介绍中提到，如果$X^{(1)},\cdots,X^{(k)}$是独立的，那么$\Sigma_{ij}=O$对任何$i\ne j$都成立，反之也成立，因此在正态总体下，假设检验就变成了以下的形式：

$H_0:\forall i\ne j,\Sigma_{ij}=O\Leftrightarrow H_1:\exist i\ne j,\Sigma_{ij}\ne O.$

由于样本均值、样本离差阵是对总体均值、自协方差矩阵的估计，因此我们也可以对样本均值和样本离差阵作同型分解。如果$H_0$成立，则$X_{(\alpha)}^{(t)}\sim N_p(\mu^{(t)},\Sigma_{tt})$且相互独立，那么似然函数就是

$L(\mu,\Sigma)=\prod_{t=1}^kL_t(\mu^{(t)},\Sigma_{tt}),$

取最大值的情况显然是$\mu^{(t)}=\bar X^{(t)},\Sigma_{tt}=A_{tt}/n$，所以似然比统计量的分子是

$\begin{aligned}
&\prod_{t=1}^n(2\pi)^{-np_t/2}|A_{tt}/n|^{-n/2}\exp\left\{-\frac12\sum_{\alpha=1}^n(X_{(\alpha)}^{(t)}-\bar X^{(t)})'\left(\frac{A_{tt}}{n} \right)^{-1}(X_{(\alpha)}^{(t)}-\bar X^{(t)}) \right\}\\
=&(2\pi)^{-np/2}\exp\left\{-\frac12\sum_{\alpha=1}^n(X_{(\alpha)}-\bar X)'\left(\frac{A}{n} \right)^{-1}(X_{(\alpha)}-\bar X) \right\}\prod_{t=1}^k\left|\frac{A_{tt}}n{} \right|^{-n/2}.
\end{aligned}$

这里的转换可以用之前常用的迹变换得出。观察分子与分母，发现其大部分是相同的，所以得到似然比统计量为

$\lambda =\frac{\prod_{t=1}^k|A_{tt}/n|^{-n/2}}{|A/n|^{-n/2}}=\left(\frac{|A|}{\prod_{t=1}^n|A_{tt}|} \right)^{n/2}\stackrel {\rm def}=V^{n/2}.$

所以我们取检验统计量为

$V=\frac{|A|}{\prod_{i=1}^k|A_{tt}|}.$

并且有结论保证，在$H_0$成立的条件下，$-b\ln V\stackrel {H_0}\to \chi^2(f)$，这里

$b=n-\frac32-\frac{p^3-\sum_{t=1}^k p_t^3}{3(p^2-\sum_{t=1}^k p_t^2)}, \\
f=\frac 12\left[p(p+1)-\sum_{t=1}^k p_t(p_t+1) \right].$

事实上$-b\ln V$是$-2\ln \lambda$的近似，故$b$也是$n$的近似，而$f$就是两个参数空间的维度之差。
2.一元数据正态性检验
回顾我们之前提到的假设检验，包括均值向量、自协方差矩阵、独立性的检验，都基于一个前提——总体是多维正态分布，如果这个正态性不满足，与三大分布相关的统计量转化、似然比统计量的表现形式都将不同于此形式，从而无法应用已有的结论。因此，本节探讨样本的正态性检验，概括起来就是，给定$n$个$p$维样本$X_{(\alpha)}$，判断总体$X$是否服从$N_p(\mu,\Sigma)$分布。
多元数据的正态性检验问题，常常转化为多个一元或二元数据的正态性检验，或者先求$X$的分量的线性组合再化为一元数据的正态性检验等。虽然我们知道，边缘分布的正态性不能推出总体分布的正态性，但是在实际应用中，这种情况并不常见，所以我们可以先将目光放在一元数据的正态性检验。
常用于一元数据检验的方法有Pearson $\chi^2$检验法（比较适合离散情形）、Kolmogorov检验法（比较适合连续情形），不过在Kolmogorov检验中我们需要得知总体的参数，即均值和方差，在实际应用中这个条件很难满足，所以我们会使用总体均值和总体方差代替，这就是Lilliefors检验。
还有一些仅适用于正态分布的检验法：偏度峰度检验法，Q-Q图和P-P图检验法、Anderson-Darling统计量检验法、Cramer-von Mises统计量检验法等。
偏度峰度法指的是，计算样本偏度和样本峰度：

$G_1=\frac{\sum(X_i-\bar X)^3}{[\sum(X_i-\bar X)^2]^{3/2}},\quad G_2=\frac{\sum(X_i-\bar X)^4}{[\sum(X_i-\bar X)^2]^2},$

在正态性成立时，近似有

$G_1\sim N\left(0,\frac{6(n-2)}{(n+1)(n+3)} \right), \\
G_2\sim N\left(3-\frac6{n-1},\frac{24n(n-2)(n-3)}{(n+1)^2(n+3)(n+5)} \right).$

很容易用Z检验找到其拒绝域。
Q-Q(Quantile Quantile)图检验法是一种图示检验法，绘制$(q_i,x_{(i)}^*)$散点图，这里$q_i=\Phi^{-1}(p_i)$是样本的$p_i$分位数，$x_{(i)}^*$是样本的$p_i$分位数，如果$X$是一元正态总体，则这些散点应该散布在一条直线上。P-P图检验法也是图示检验，绘制的数据点是$(p_i,F(x_{(i)}^*))$，其中$p_i$是经验分布函数$F_n(x)$在$x_{(i)}^*$上的值，$F(x_{(i)}^*)$是$\Phi(x)$在$x_{(i)}^*$上的值。在实际应用Q-Q图检验和P-P图检验时，$x_{(i)}^*$要先选好。
Anderson-Darling $A^2$检验（AD检验）的检验统计量是

$A^2=n\int_{-\infty}^\infty \frac{(F_n(x)-\Phi(x))^2}{\Phi(x)(1-\Phi(x))}{\rm d}\Phi(x),$

这里$[\Phi(x)(1-\Phi(x))]^{-1}$是权重函数，如果权重函数取$1$，就得到Cramer-von Mises $W^2$检验的检验统计量

$W^2=n\int_{-\infty}^\infty (F_n(x)-\Phi(x))^2{\rm d}\Phi(x).$

结合Kolmogorov-Smirnov统计量$D=\sup|F_n(x)-\Phi(x)|$，这三个统计量都是原假设成立时不能过大的，依赖于一个概率表值来检验原假设是否应该被接受。不过，这三种检验方式适用于各种假设检验，只要将表达式中的$\Phi(x)$换成对应的分布函数即可。
3.多元数据的正态性检验
对于二元数据，存在一种粗糙的检验方法：等概椭圆检验法。其理论基础是二维随机向量$X$如果来自于正态总体，则其概率密度函数等高线应该是一个椭圆，即$X\sim N_2(\mu,\Sigma)$时，应有

$f(x_1,x_2)=a\Leftrightarrow (X-\mu)'\Sigma^{-1}(X-\mu)=b^2.$

所以我们计算二元数据$X_{(i)}$到$\bar X$的马氏距离$D_i=(X_{(i)}-\bar X)'S^{-1}(X_{(i)}-\bar X)$，在给定数值$p_0$下，$D_i\le p_0$的频率应该和某一个定值比较接近，这个定值可以通过查表获得。由于这是一种比较粗糙的方法，我们在实际应用中会使用更为正式的方法。
现在介绍$p$维数据$\chi^2$统计量的Q-Q图检验法，我们将假设确定为参数已知的，即

$H_0:X\sim N_p(\mu,\Sigma)\Leftrightarrow H_1:X\nsim N_p(\mu,\Sigma).$

由于在正态性假设$H_0$成立的前提下，样本$X$到中心$\mu$的马氏距离存在以下关系：

$D^2=(X-\mu)'\Sigma^{-1}(X-\mu)\sim \chi^2(p),$

所以我们可以直观地想到验证样本的马氏距离是否具有这样的关系。因此，我们计算样本$X_{(\alpha)}$到$\mu$的马氏距离$D_{\alpha}^2=(X_{(\alpha)}-\mu)'\Sigma^{-1}(X_{(\alpha)}-\mu)$，并对$D_\alpha^2$进行排序得到次序统计量$D_{(\alpha)}$，计算其经验分布函数，这样有了经验分布函数与$\chi^2(p)$分布的分布函数后，就可以绘制Q-Q图或者P-P图。
在实际应用中，我们往往不知道$\mu,\Sigma$的值，所以会用样本均值$\bar X$和样本协方差阵$A/(n-1)$代替，得到的Q-Q图或P-P图应该是一条通过原点、斜率为1的直线，如果是这样，就可以接受正态性假设，否则应当拒绝。
回顾总结


正态总体的独立性检验，我们一般会取检验统计量为

$V=\frac{|A|}{\prod_{t=1}^k |A_{tt}|}.$

当$n\to \infty$时，有$-b\ln V\to \chi^2(f)$，这里

$b=n-\frac32-\frac{p^3-\sum_{t=1}^k p_t^3}{3(p^2-\sum_{t=1}^k p_t^2)},\\
f=\frac{p(p+1)}{2}-\sum_{t=1}^k\frac{p_k(p_k+1)}{2}.$


一元总体的正态性检验有很多方法，如K-S检验、A-D检验、Cramer-von Mises检验，但K-S检验的效果一般，A-D检验的效果比较好，其检验统计量是

$A^2=n\int_{-\infty}^\infty \frac{(F_n(x)-\Phi(x))^2}{\Phi(x)(1-\Phi(x))}{\rm d}\Phi(x).$


Q-Q图是分位数图，首先选定一组分位数间隙$x_{(i)}^*$，然后在样本中寻找相应分位数，在总体中也寻找相应分位数，将分位数绘制成散点图，观察其是否位于一条直线上。


P-P图是累计分布图，首先选定一组分位数间隙$x_{(i)}^*$，然后绘制经验分布函数与总体分布函数在$x_{(i)}^*$处的取值，将两个取值绘制成散点图，观察其是否位于一条直线上。


多元总体的正态性检验采用$\chi^2$统计量的Q-Q图检验法，计算样本到中心$\bar X$的马氏距离并排序，用Q-Q图判断是否属于$\chi^2(p)$分布，或用K-M检验法。马氏距离的定义如下：

$D_\alpha=(X_{(\alpha)}-\bar X)'S^{-1}(X_{(\alpha)}-\bar X).$

课堂笔记

多元正态分布的密度函数

每一个分量都服从正态分布
分量的线性组合仍然是正态分布
如果协方差矩阵是对角阵，则分量是相互独立的、服从正态分布的随机变量（对于正态分布而言，不线性相关等价于独立）。

正态性的识别
图方法
Q-Q图：点是否在直线的附近

x 理论分位数
y 样本分位数：修正的计算公式 $\frac{i-0.357}{n+0.25}$
结论：$X_{(i)}=\mu+\sigma Z_{(i)}$


$\alpha$分位数：$P(x<=X_{\alpha}) =\alpha$

R语言实现：qqnorm() qqline()

正态W检验：shapiro.test()

原假设：正态成立

W越靠近1，越好

概率小，拒绝正态性假设

对数变换增强正态性：

Cox-Box变换（考试会考）

一、业务场景：

1. 一个汽车销售公司，其客户来店消费金额是否符合正态分布？

答：这个问题可以抽象为统计学的统计推断中的假设检验部分的正态性检验。

2. 如何模拟这些数据的函数特征，怎么看拟合的好不好？

答：这是个拟合问题，视情况用线性拟合和多项式拟合来拟合。通过拟合打分看拟合效果。

3. 这个具体函数能否给出来？

答：可以。


二、下面分四部分来用代码解决上述问题

1. 对数据做正态性判断

2. 对数据做多元线性回归

3. 对数据做多元多项式回归

4. 输出线性回归和多项式回归对应的具体函数（函数包括：组合特征项，特征系数和常数项）


1. 正态性检验（Normality Test）

#三种方法：

#1. 生成数据为浮点数，所以用密度，如果数据是单词出现频率等整数则用直方图比较，根据图形肉眼比较。

#2. 用shapiro-wilk norm test，p-value的值大于0.05可视为满足正态性。

#3. 用qqplot做正态性检验，解读方法是如果qqplot所划散点和直线越接近则是正态分布，否则不是。


#代码

## Generate two data sets

## First Normal, second from a t-distribution

words1 = rnorm(100); words2 = rt(100, df=3)


## Have a look at the densities

plot(density(words1));plot(density(words2))


## Perform the shapiro-wilk test

shapiro.test(words1); shapiro.test(words2)


## Plot using a qqplot

qqnorm(words1);qqline(words1, col = 2)

qqnorm(words2);qqline(words2, col = 2)


参考：

https://stats.stackexchange.com/questions/3136/how-to-perform-a-test-using-r-to-see-if-data-follows-normal-distribution

https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.stats.shapiro.html



2. 用scikit-learn做多元线性回归

# R-squared表示模型的拟合度好坏，取值为0-1，越接近1越好。

from sklearn.linear_model import LinearRegression

X = [[6, 2], [8, 1], [10, 0], [14, 2], [18, 0]]

y = [[7], [9], [13], [17.5], [18]]

model = LinearRegression()

model.fit(X, y)

X_test = [[8, 2], [9, 0], [11, 2], [16, 2], [12, 0]]

y_test = [[11], [8.5], [15], [18], [11]]

predictions = model.predict(X_test)

for i, prediction in enumerate(predictions):

    print('Predicted: %s, Target: %s' % (prediction, y_test[i]))


# 或者用summary(r)查看R-squared，分别查看训练集和测试集

print('R-squared: %.2f' % model.score(X_test, y_test))

print('R-squared: %.2f' % model.score(X, y))


参考：http://blog.csdn.net/SA14023053/article/details/51703204


3. 对数据做多元多项式回归

import numpy as np

import matplotlib.pyplot as plt


from sklearn.linear_model import Ridge

from sklearn.preprocessing import PolynomialFeatures

from sklearn.pipeline import make_pipeline


def f(x):

    """ function to approximate by polynomial interpolation"""

    return x * np.sin(x)


# generate points used to plot

x_plot = np.linspace(0, 10, 100)


# generate points and keep a subset of them

x = np.linspace(0, 10, 100)

rng = np.random.RandomState(0)

rng.shuffle(x)

x = np.sort(x[:20])

y = f(x)


# create matrix versions of these arrays

X = x[:, np.newaxis]

X_plot = x_plot[:, np.newaxis]


colors = ['teal', 'yellowgreen', 'gold']

lw = 2

plt.plot(x_plot, f(x_plot), color='cornflowerblue', linewidth=lw,

         label="ground truth")

plt.scatter(x, y, color='navy', s=30, marker='o', label="training points")



for count, degree in enumerate([3, 4, 5]):

    model = make_pipeline(PolynomialFeatures(degree), Ridge())

    model.fit(X, y)

    y_plot = model.predict(X_plot)
#拟合出来的是函数是个靠近值不是经过值，所以可得到训练集的拟合程度。

    print model.named_steps['ridge'].coef_

    print('score: %s' % (model.score(X, y)))

    plt.plot(x_plot, y_plot, color=colors[count], linewidth=lw,

             label="degree %d" % degree)


plt.legend(loc='lower left')

plt.show()


参考：http://scikit-learn.org/stable/auto_examples/linear_model/plot_polynomial_interpolation.html


4. 输出线性回归和多项式回归对应的具体函数

具体函数：函数的组成包括组合特征项，特征系数和常数项。类似y = const + ax1 + bx2 + cx1^2 + dx1x2 + ex2^2。

对于一元/多元，线性/多项式等不同情况，我们只要把函数各组成部分找出来，就能把具体函数列出来。

特征项：

1> 对于单个输入特征x1，线性，度为1, 拟合的函数为y = 1+x1

2> 对于单个输入特征x1，度为2，拟合的函数为y = (1+x1)^2 = 1 + 2x1 + x1^2

3> 对于n个(以2为例)输入特征x1和x2，度为2，拟合的函数为 y = (1+x1+x2)^2 = 1 + 2(x1+x2) + (x1+x2)^2 = 1 + 2x1 + 2x2 + x1^2 + 2x1x2 + x2^2

去掉系数就是组合特征项。

特征系数：model.coef_可得到所有组合特征项的系数。

常数项：用model.predict(0, 0)得到，也就是y - ax1 - bx2 = 常数项


参考：http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html