【数据建模学习实录】——数据处理与模型拟合

df = pd.DataFrame({'From_To': ['LoNDon_paris', 'MAdrid_miLAN', 'londON_StockhOlm',
                               'Budapest_PaRis', 'Brussels_londOn'],
              'FlightNumber': [10045, np.nan, 10065, np.nan, 10085],
              'RecentDelays': [[23, 47], [], [24, 43, 87], [13], [67, 32]],
                   'Airline': ['KLM(!)', '<Air France> (12)', '(British Airways. )',
                               '12. Air France', '"Swiss Air"']})
df

	From_To	FlightNumber	RecentDelays	Airline
0	LoNDon_paris	10045.0	[23, 47]	KLM(!)
1	MAdrid_miLAN	NaN	[]	<Air France> (12)
2	londON_StockhOlm	10065.0	[24, 43, 87]	(British Airways. )
3	Budapest_PaRis	NaN	[13]	12. Air France
4	Brussels_londOn	10085.0	[67, 32]	"Swiss Air"数据的

（3）插值模型

1. 线性插值

import numpy as np
#数据准备
X = np.arange(-2*np.pi, 2*np.pi, 1) # 定义样本点X，从-pi到pi每次间隔1
Y = np.sin(X) # 定义样本点Y，形成sin函数
new_x = np.arange(-2*np.pi, 2*np.pi, 0.1) # 定义插值点
# 进行样条插值
import scipy.interpolate as spi
# 进行一阶样条插值
ipo1 = spi.splrep(X, Y, k=1)  # 样本点导入，生成参数
iy1 = spi.splev(new_x, ipo1)  # 根据观测点和样条参数，生成插值

2. 三次样条插值

# 进行三次样条拟合
ipo3 = spi.splrep(X, Y, k=3)  # 样本点导入，生成参数
iy3 = spi.splev(new_x, ipo3)  # 根据观测点和样条参数，生成插值

3.拉格朗日插值

def lagrange(x0,y0,x):
    y=[]
    for k in range(len(x)):
        s=0
        for i in range(len(y0)):
           t=y0[i]
           for j in range(len(y0)):
              if i!=j:
                t*=(x[k]-x0[j])/(x0[i]-x0[j])
           s+=t
        y.append(s)
    return y

2.3 数据的规约

两种典型的规约方式

（1） min-max规约：

$eq?%7Bx_%7Bnew%7D%7D%20%3D%20%5Cfrac%7B%7Bx%20-%20%5Cmin%20%28x%29%7D%7D%7B%7B%5Cmax%20%28x%29%20-%20%5Cmin%20%28x%29%7D%7D$

所有的属性都会被规约到【0，1】的范围内，但是对异常值表现明显，即如果出现异常值，这个数据的分布就是有偏的。

（2）Z-score规约：

$eq?%7Bx_%7Bnew%7D%7D%20%3D%20%5Cfrac%7B%7Bx%20-%20%5Cbar%20x%7D%7D%7B%7Bstd%28x%29%7D%7D$

本质上，一列数据减去其均值再除以标准差，如果这一列数据近似服从正态分布，这个过程就是化为标准正态分布的过程。Z-score规约和min-max规约往往不是二者取其一，有时候两个可以组合起来用。

3. 常见的统计分析模型

3.1 回归分析与分类分析

回归分析中因变量是连续变量，如工资、销售额；而分类分析中因变量是属性变量，如判断邮件“是or否”为垃圾邮件。

（1）回归分析：在学生刚入大学时，大学通常会举办一个大学入学考试，该成绩虽然没什么用，但是能让刚入学的学生保持持续学习的警示。现在，我们想分析大学成绩GPA是否与入学成绩以及高考成绩有关？

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
from IPython.display import display
import statsmodels.api as sm

# 加载数据
gpa1=pd.read_stata('./data/gpa1.dta')

# 在数据集中提取自变量（ACT为入学考试成绩、hsGPA为高考成绩）
X1=gpa1.ACT   # x1 是只有 一个入学成绩
X2=gpa1[['ACT','hsGPA']]  # x2是含有入学成绩和高考成绩
# 提取因变量
y=gpa1.colGPA

# 为自变量增添截距项
X1=sm.add_constant(X1)
X2=sm.add_constant(X2)
display(X2)


# 拟合两个模型
gpa_lm1=sm.OLS(y,X1).fit()
gpa_lm2=sm.OLS(y,X2).fit()

# 输出两个模型的系数与对应p值
p1=pd.DataFrame(gpa_lm1.pvalues,columns=['pvalue'])
c1=pd.DataFrame(gpa_lm1.params,columns=['params'])


p2=pd.DataFrame(gpa_lm2.pvalues,columns=['pvalue'])
c2=pd.DataFrame(gpa_lm2.params,columns=['params'])
display(c1.join(p1,how='right'))
display(c2.join(p2,how='right'))

注意：X1=gpa1.ACT # x1 是只有一个入学成绩
X2=gpa1[['ACT','hsGPA']] # x2是含有入学成绩和高考成绩

我们分析结果首先看第一个拟合模型结果， 1.3899927e-02 小于 0.05 ，因此入学成绩能显著影响大学成绩

但从模型II可以看到，高考成绩haGPA的pvalue为0.000005远小于0.05，说明高考成绩haGPA显著影响大学成绩GPA，但大学入学考试成绩ACT的pvalue为0.38大于0.05，不能说明大学入学考试成绩ACT显著影响大学成绩GPA。

显然，模型1 和模型2 的结果是相悖的，这是为什么呢？

对于模型1和2，两者的区别就是，影响因素的选择不同，模型1只有一个入学成绩，而模型2 含有入学成绩和高考成绩

所以出现相悖结果的原因就是高考成绩和大学入学考试成绩的相关性也很高，但是高考考试成绩对GPA的影响更加显著，就显得大学成绩的影响不那么显著了。

3.2 假设检验

这里需要注意的是，在实际应用中，我们应该考虑到数据的复杂性，考虑同时使用多种方法进行判断

# 接下来，我们在python中定义一个函数，将概率图、Shapiro-Wilk test、D'Agostino's K-squared test结合在一起。
data_small = stats.norm.rvs(0, 1, size=30) # 小样本正态性数据集
data_large = stats.norm.rvs(0, 1, size=6000) # 大样本正态性数据集
# 定义一个正态性检验函数，它可以输出：
## 正态概率图
## 小样本Shapiro-Wilk检验的p值
## 大样本D'Agostino's K-squared检验的p值

from statsmodels.stats.diagnostic import lilliefors
from typing import List

def check_normality(data: np.ndarray, show_flag: bool=True) -> List[float]:
    """
    输入参数
    ----------
    data : numpy数组或者pandas.Series
    show_flag : 是否显示概率图
    Returns
    -------
    两种检验的p值；概率图
    """

    if show_flag:
        _ = stats.probplot(data, plot=plt)
        plt.show()

    pVals = pd.Series(dtype='float64')
    # D'Agostino's K-squared test
    _, pVals['Omnibus'] = stats.normaltest(data) 

    # Shapiro-Wilk test
    _, pVals['Shapiro-Wilk'] = stats.shapiro(data)

    print(f'数据量为{len(data)}的数据集正态性假设检验的结果 : ----------------')
    print(pVals)
check_normality(data_small,show_flag=True)

数据量为30的数据集正态性假设检验的结果 : ----------------

Omnibus 0.703004

Shapiro-Wilk 0.898077

数据量为6000的数据集正态性假设检验的结果 : ----------------

Omnibus 0.980014

Shapiro-Wilk 0.992427