Datawhale 数学建模导论第六章学习笔记

数据处理与拟合模型

数据与大数据

数据

• 数据的定义：具有描述与传递信息功能的符号记录。
• 信息的模态定义：信息的不同的形式。

       例如：

• 数值类数据，例如结构化的excel表格和SQL文件。
• 文本类数据，例如新闻报道、微博评论、餐饮点评等文字。
• 图像类数据，以一定尺寸的黑白或彩色图像在计算机内存储。
• 音频类数据，例如音乐、电话录音等。
• 信号类数据，例如地震波的波形、电磁波信号、脑电信号等。

不同形式的信息均可数化 ，而不同模态的数据又往往能联合起来对同一事物进行描述。在此基础上建立的模型称之为多模态模型。

以体量区分数据：

• 大数据往往是以TB甚至PB、ZB来衡量。
• 小数据是在G以内的数据。

从数据的体量的角度来看（处理数据时建立模型的根本衡量角度），建立模型不能随意地“大摆身手”， 面对不同的数据，能够使用最合适的方法最为重要。

数据科学的研究对象

 

大数据科学研究的不仅仅是数据分析，它还包括了：

• 数据的获取和存储：包括爬虫、软件定义存储、硬件存储有关背景知识等。
• 数据的处理：包括分布式计算、并行计算、数据流等知识，以及Hadoop、Spark等大数据框架。
• 数据的分析：包括统计学、数据挖掘与机器学习、计算机视觉、自然语言处理等内容，重在挖掘数据中的模式与知识。
• 数据的管理：现代数据库系统及其架构等内容。
• 数据的应用：数据可视化、数据相关软件的开发、报表分析以及如何将数据挖掘得到的结果还原为实际问题的解决方案。

 数据的预处理

数据预处理的原因

在现实问题中得到的原始数据常常混乱或者是不全，数据和特征决定了效果的上限，而模型只是逼近这个上限而已，可见数据预处理是十分重要的。

 表格数据的稀疏划分

表格数据的列表示这个表格数据的一个属性，属性又可以称之为维（度）。有几列就是有几个属性，就是有几个维度。

属性分为：连续属性，离散属性。例如在图中，“品牌类型”这个属性只有{1,2,3}三个有限的取值，我们称之为离散。当然，这个取值也可以不是数字，比如{汽车，火车，飞机}等。

稀疏：列数超过行数的1/2，是行数的3倍就是严重稀疏。

稀疏另一种定义：表格里面很多都是空白的或者全是0。

数据有的问题包括：空缺、重复和一些异常记录。解决措施：

• 重复：直接将其删除即可。
• 空缺：通过观察缺失率处理问题

1. 如果存在缺失的数据项（数表一行称作一个数据项）占比较少（大概5%以内）这个时候如果问题允许可以把行删掉。
2. 如果缺失率稍微高一点（5%-20%）左右就可以使用填充、插值的方法去处理，有关插值的方法会在下一节探讨，而填充方法包括常数填充、均值填充等方法；
3. 如果缺失率还高一些（20%-40%）那么就需要用预测方法例如机器学习去填充缺失数据了；但如果一行数据有50%以上都是缺失的，如果条件允许，我们可以把这一列都删掉（当然凡事都有例外，见机行事）。

pandas在数据预处理的应用例子：

'''（1）创建pandas dataframe'''
import pandas as pd
import numpy as np
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"']})

print(df.describe())
'''（2）FlightNumber列中有某些缺失值，缺失值常用nan表示，请在该列中添加10055与10075填充该缺失值。'''
df['FlightNumber'] = df['FlightNumber'].interpolate().astype(int)
print(df.describe())
'''（3）由于列From_To 代表从地点A到地点B，因此可以将这列拆分成两列，并赋予为列From与To。'''
temp = df['From_To'].str.split("_", expand=True)
temp.columns = ['From', 'To']

'''（4）将列From和To转化成只有首字母大写的形式。'''
temp['From'] = temp['From'].str.capitalize()
temp['To'] = temp['To'].str.capitalize()

'''（5）将列From_To从df中去除，并把列From和To添加到df中'''
df.drop('From_To', axis=1, inplace=True)
df[['From', 'To']] = temp
print(df[['From', 'To']] )

'''（6）清除列中的特殊字符，只留下航空公司的名字。'''
df['Airline'] = df['Airline'].str.extract(r'([a-zA-Z\s]+)', expand=False).str.strip()
print(df[['Airline']] )
'''（7）在 RecentDelays 列中，值已作为列表输入到 DataFrame 中。
我们希望每个第一个值在它自己的列中，每个第二个值在它自己的列中，依此类推。
如果没有第 N 个值，则该值应为 NaN。将 Series 列表展开为名为 的 DataFrame delays，
重命名列delay_1，delay_2等等，并将不需要的 RecentDelays 列替换df为delays。'''
delays = df['RecentDelays'].apply(pd.Series)
delays.columns = ['delay_%s' % i for i in range(1, len(delays.columns)+1)]
df = df.drop('RecentDelays', axis=1).join(delays, how='left')
print(df)
'''（8）将delay_i列的控制nan都填为自身的平均值。'''
for i in range(1, 4):
    df[f'delay_{i}'] = df[f'delay_{i}'].fillna(np.mean(df[f'delay_{i}']))

'''（9）在df中增加一行，值与FlightNumber=10085的行保持一致'''
df = df.append(df.loc[df['FlightNumber'] == 10085, :], ignore_index=True)

'''（10）对df进行去重，由于df添加了一行的值与FlightNumber=10085的行一样的行，因此去重时需要去掉。'''
df = df.drop_duplicates()
print(df)

 常见统计分析模型

回归分析和分类分析

• 共同点：是基于统计模型的统计分析方法，都研究因变量和自变量之间存在的潜在关系，并通过统计模型的形式将这些潜在关系显式的显示出来。
• 不同点：回归分析：因变量是连续变量，如工资、销售额；分类分析：因变量是属性变量，如判断邮件“是or否”为垃圾邮件。

到底是采用回归分析还是分类分析的判断依据就是因变量是连续变量还是属性变量。 

假设检验 

假设检验：根据样本信息与已知信息，对一个描述总体性质的命题进行“是或否”的检验与回答。假设检验验证的不是样本本身的性质，而是样本所在总体的性质。

假设检验在建模过程的用处：

1. 在建模中对数据进行探索性的信息挖掘，提供了选择模型的依据。
2. 在建模完成后，可以检验该模型的有效性。

分类：

1. 参数假设检验：假设是关于总体的一个参数或是参数的集合；
2. 非参数假设检验：假设不能用一个参数集合表示，例如：正态性检验。 

• 正态性检验：如果数据是正态分布的， 我们应该采用正态性检验。
• 判断数据的正态性方法：

1.  可视化判断-正态分布概率图
2. Shapiro-Wilk检验
3. D'Agostino's K-squared检验

实际处理数据的过程中，我们需要综合使用多种方法进行判断数据的正态性。如果不同方法得出的结论不同，我们就要观察数据的特征，找到结论不一致的原因。如：若Shapiro-Wilk test显著（非正态），D'Agostino's K-squared test不显著（正态），则有可能是因为样本量较大，或者样本中存在重复值现象，如果事实确实如此，那么我们就应该采纳D'Agostino's K-squared test的结论而非Shapiro-Wilk test的结论。（在这里，我们假设数据服从正态分布，p<0.05则拒绝假设，换个说法：p<0.05数据不服从正态分布，如果p>=0.05说明没有证据说明数据不服从正态分布）。

举例的对应代码如下：

# 接下来，我们在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
dtype: float64

check_normality(data_large,show_flag=False) # 当样本量大于5000，会出现警告
数据量为6000的数据集正态性假设检验的结果 : ----------------
Omnibus         0.980014
Shapiro-Wilk    0.992427
dtype: float64

e:\anaconda\envs\ml\lib\site-packages\scipy\stats\morestats.py:1760: UserWarning: p-value may not be accurate for N > 5000.
  warnings.warn("p-value may not be accurate for N > 5000.")

•  单组样本均值假定检验：样本所代表的总体的xxx是否能与总体的xxx有显著区别。例如：

  必胜中学里，陈老师班结束了一次英语考试。由于班级人数较多，短时间内很难完成批改与统计，陈老师又很想知道此次班级平均分与级长定的班级均分137的目标是否有显著区别，于是他随机抽取了已经改好的10名同学的英语成绩：

  136,136,134,136,131,133,142,145,137,140

  问：陈老师可以认为此次班级平均分与级长定的班级均分137的目标没有显著区别吗？

  很明显，这就是一个典型的单组样本均值假定的检验，比较的是这个样本（10个同学的英语成绩）所代表的总体均值（班级英语成绩均值）是否与参考值137相等。

•   单组样本均值假定检验的分类：

1. 单样本t检验
2. wilcoxon检验 

• 两组样本的均值相等性检验：

独立：两个样本中，一个样本中的受试不能影响另一个样本中的受试。

但是实际中，事件是否“独立”是难以判断的，所以一般只能在抽样的时候使用更加科学的抽样方法，尽可能地避免样本间的相互影响。换句话说，虽然以是否“独立”作为选择检验方法的依据不靠谱，但是依据可以将其作为初步判断的依据，所以对于两组样本的均值相等性检验，我们通常使用双样本t检验，Mannwhitneyu秩和检验，成对检验。

• 双样本t检验：两个样本的总体都服从正态分布。
• Mannwhitneyu秩和检验：两个样本的总体都不服从正态分布。

成对：两个样本分别为同一个受试个体不同时间的受试结果，这两个样本是“成对”的，是彼此紧密相连的。而对这样两个样本进行均值比较检验，就是成对检验。

• 成对检验分为：

1. 成对t检验：若总体服从正态分布
2. 成对wilcoxon秩和检验：若总体不服从正态分布。

综上所述：先判断是独立还是成对，然后再判断数据是否为正态分布。 

• 方差分析-多组样本间的均值相等性检验：

 既然我们检验的是不同总体的均值是否相等，那么观察各样本的样本均值的“差异程度”一定是非常自然且合理的想法，如果各样本的均值差异很大，那么它们的总体均值也有很大可能存在差异。但是这并不够，样本的样本内差异程度也十分重要。样本内差异程度越大，“偶然性”越大，我们越难以判断两个不相等的均值是否真的不相等。因此，我们需要综合这两个评判指标。最简单的方法就是两者相除，即样本间均值的“差异程度”除以样本内差异程度。这就是方差分析最根本的思想。

举个例子：小红的考试均分是91，小刚的考试均分是89。我们假设一个非常极端的情况：他们的标准差都是0，即小红每次考试都是91，小刚每次考试都是89，那么我们似乎可以很容易地判断出，两个人的均分确实存在明显的差异；但是，如果他们的标准差都很大，高达6（方差就是36），即他们的成绩都很不稳定，这次小红考79、小刚考93，下次小红考94、小刚考70。在高达36的方差下，2分的均值差似乎没有什么说服力了。

尽管在大样本下，非正态性数据的方差分析也是稳健的，但是在小样本下，对非正态性数据做方差分析还是可能存在误差。此时，我们可以使用kruskalwallis检验。

综上所述，建模前需要对数据进行处理：

1. 数据预处理：对数据进行填补空缺，去重，使得数据变得更可用。
2. 假设检验：提供选择模型的依据。

随机过程与随机模拟

•  随机过程：描述一段时间内，随机变量的分布的变化。随机过程就是在随机变量的基础上加入了时间维度（值得注意的是，时间维度不是随机变量，只是普通变量）。
• 案例分析：

 

为了改善道路的路面情况（道路经常维修，坑坑洼洼），因此想统计一天中有多少车辆经过，因为每天的车辆数都是随机的，一般来说有两种技术解决这个问题：

（1） 在道路附近安装一个计数器或安排一个技术人员，在一段长时间的天数（如365天）每天24h统计通过道路的车辆数。

（2） 使用仿真技术大致模拟下道路口的场景，得出一个近似可用的仿真统计指标。

由于方案（1）需要花费大量的人力物力以及需要花费大量的调研时间，虽然能得出准确的结果，但是有时候在工程应用中并不允许。因此，我们选择方案（2），我们通过一周的简单调查，得到每天的每个小时平均车辆数：[30, 20, 10, 6, 8, 20, 40, 100, 250, 200, 100, 65, 100, 120, 100, 120, 200, 220, 240, 180, 150, 100, 50, 40]，通过利用平均车辆数进行仿真。

# 模拟仿真研究该道路口一天平均有多少车经过
import simpy

class Road_Crossing:
    def __init__(self, env):
        self.road_crossing_container = simpy.Container(env, capacity = 1e8, init = 0)
    
def come_across(env, road_crossing, lmd):
    while True:
        body_time = np.random.exponential(1.0/(lmd/60))  # 经过指数分布的时间后，泊松过程记录数+1
        yield env.timeout(body_time)  # 经过body_time个时间
        yield road_crossing.road_crossing_container.put(1)

hours = 24  # 一天24h
minutes = 60  # 一个小时60min
days = 3   # 模拟3天
lmd_ls = [30, 20, 10, 6, 8, 20, 40, 100, 250, 200, 100, 65, 100, 120, 100, 120, 200, 220, 240, 180, 150, 100, 50, 40]   # 每隔小时平均通过车辆数
car_sum = []  # 存储每一天的通过路口的车辆数之和
print('仿真开始：')
for day in range(days):
    day_car_sum = 0   # 记录每天的通过车辆数之和
    for hour, lmd in enumerate(lmd_ls):
        env = simpy.Environment()
        road_crossing = Road_Crossing(env)
        come_across_process = env.process(come_across(env, road_crossing, lmd))
        env.run(until = 60)  # 每次仿真60min
        if hour % 4 == 0:
            print("第"+str(day+1)+"天，第"+str(hour+1)+"时的车辆数：", road_crossing.road_crossing_container.level)
        day_car_sum += road_crossing.road_crossing_container.level
    car_sum.append(day_car_sum)
print("每天通过交通路口的的车辆数之和为：", car_sum)

仿真开始：
第1天，第1时的营业额： 80
第1天，第5时的营业额： 44
第1天，第9时的营业额： 996
第1天，第13时的营业额： 989
第1天，第17时的营业额： 582
第1天，第21时的营业额： 582
第2天，第1时的营业额： 81
第2天，第5时的营业额： 65
第2天，第9时的营业额： 833
第2天，第13时的营业额： 917
第2天，第17时的营业额： 496
第2天，第21时的营业额： 585
第3天，第1时的营业额： 85
第3天，第5时的营业额： 89
第3天，第9时的营业额： 888
第3天，第13时的营业额： 1048
第3天，第17时的营业额： 554
第3天，第21时的营业额： 734
每天商店的的营业额之和为： [11491, 11014, 11572]

数据可视化 

• 数据可视化：更加强调使用可视化的图表去表达数据中的信息，如数据分析中的结论。
• 数据可视化分类：

1. 在分析过程中的数据可视化：强调辅助分析，即图表不是很在意观感如何，只要自己能看懂，自己从图表中得到有用的结论并辅助分析的流程其实就可以了。
2. 分析结果表达中的数据可视化：强调在阅读数据分析结论的人能更好地知道分析结论是什么，因此相对于分析过程中的数据可视化来说更加美观、表达的信息更全面。

• 好的数据可视化图表具有的条件：

1. 图表展示的信息全面且无歧义
2. 图表表达的信息越多、越全面越好
3. 通俗易懂，不能太专业

Python三大数据可视化工具库：

1.  Matplotlib

2. Seaborn
3. Plotnine

 基本图标Quick Start

1. 类别型图表：类别型图表一般表现为：X类别下Y数值之间的比较，因此类别型图表往往包括：X为类别型数据、Y为数值型数据。类别型图表常常有：柱状图、横向柱状图（条形图）、堆叠柱状图、极坐标的柱状图、词云、雷达图、桑基图等等。
2. 关系型图表：关系型图表一般表现为：X数值与Y数值之间的关系，如：是否是线性关系、是否有正向相关关系等等。一般来说，关系可以分为：数值型关系、层次型关系和网络型关系。
3. 分布型图表：所谓数据的分布，其实就是数据在哪里比较密集，那里比较稀疏，描述数据的密集或者稀疏情况实际上可以用频率或者概率。

插值模型

插值方法可以用于处理缺失值，也可以进行数据填充。比如，我获取了按日的数据，需要按小时对数据进行填充，这个时候就可以用插值去做填充任务。这里介绍几种比较常见的插值方法。

常见插值方法：

1. 线性插值法
2. 三次样条插值
3. 拉格朗日插值 

《Python数学建模算法与应用 司守奎 》这本书对这三种方法有详细地描述，读者可以拜读此书。