首先，按照国际惯例，我们从最最简单的机器学习方案（线性回归）入手，这种简易的方案必须建立在一条基本的假设上，预报变量的分布基本服从正态分布，因此我选择了温度作为这次的预报量，（暂时不去考虑如降水量这类分布状况比较随机的量）。
同时，为了使得逐小时的温度变化具有一定的统计特征，下面使用到的观测与模式预报量都进行了如下处理：
1.获取“广东省”区域内的同一时次目标变量的所有水平空间点集。
2.将点集求平均，代表目标变量目标时次的值。
最后说明一下本章用到的资料：
1.2017年4月份逐小时的地面站观测资料（文本格式）
2.2017年4月份华南mars3km数值模式预报产品，每天00时刻（UTC）的预报产品，每次00时只取00～48小时预报时效的产品文件。（二进制顺序存取文件）
3.由于气象资料的数据太庞大，本章不提供原始的数据资料，不过我提供了处理后的应用于本章的资料集，具体操作参见后文提到的泡菜。
代码汇总于序章中提到的github中
本章的大致结构来源于我网上找到的一篇教程，的第一和第二章，有兴趣的链接在这里，链接
一. 对OBSReader模块的说明
首先，原始的观测数据资料是这样的方式命名的：
image.png
其中SYNOP文件类型，后面8位数字是站点观测对应的时间信息，如17年04月01日00时（17040100）。
里面的结构比较凌乱，既有表头，各类变量名称又是按照不同站点不同层次存放的，而且一般是全国范围的里面有很多本实验不关注区域的信息类似这样：
image.png
于是我写了一段批处理的python脚本将我们需要的单层观测站点信息提取成这样的文件：
image.png
标准的pandas快速识别格式，并只包含广东附近的站点，文本体积大大减小。
具体的处理脚本在这里（！！！！）
然后，通过OBSReader模块将指定的文本数据提取出来，具体实现如下：
# coding: utf-8
import pandas as pd #pandas主模块
import dateparser #时间类生成模块
import datetime #时间类型模块
import os #操作系统模块
import re #正则表达式模块
import geopandas as gpd #地理信息模块
from shapely.geometry import Point #用于创建包含地理信息的点阵
以上是一些基本库载入，不多赘述了。
OBSTIME = '(\d+)(\.\w+)'
P = re.compile("(SYNOP)%s" %(OBSTIME))
GUANGDONG = gpd.read_file('/home/wenqs/project_mlwf/read_binary/data/XZQ_D.shp')
class obs_from_file(object):
def __init__(self,filename,nav=-888888.0,sep=' ',incr=1):
self.filename = filename
self._find_filegroup(incr)
self._read_data(nav,sep)
self._compute_var()
def _find_filegroup(self,incr):
path = os.path.dirname(self.filename)
file = os.path.basename(self.filename)
m = P.search(file)
tail = m.group(3)
flag = m.group(1)
otime=m.group(2)
otime = dateparser.parse('20%s/%s/%s-%s:00:00' % (otime[:2],otime[2:4],otime[4:6],otime[6:8]))
otime_b = otime - datetime.timedelta(hours=incr)
otime_a = otime + datetime.timedelta(hours=incr)
self.filename_b = path+os.sep+flag+otime_b.strftime('%y%m%d%H')+tail
self.filename_a = path+os.sep+flag+otime_a.strftime('%y%m%d%H')+tail
def _read_data(self,nav,sep):
data = pd.read_csv(self.filename,sep=sep,na_values=nav,index_col='station_id')
if os.path.isfile(self.filename_b):
data_b = pd.read_csv(self.filename_b,sep=sep,na_values=nav,index_col='station_id')
else:
data_b = data.copy()
if os.path.isfile(self.filename_a):
data_a = pd.read_csv(self.filename_a,sep=sep,na_values=nav,index_col='station_id')
else:
data_a = data.copy()
data_a.columns=data_a.columns.map(lambda x : x[:]+'_a')
data_b.columns=data_b.columns.map(lambda x : x[:]+'_b')
df = pd.concat([data,data_b,data_a],axis=1)
df=df[:][pd.notnull(df['longitude'])]
need_to_drop=[]
for col in list(df):
if col.endswith('_a') or col.endswith('_b'):
need_to_drop.append(col)
continue
df_ab = df.copy()
missing_vals = pd.isnull(df[col])
df[col][missing_vals] = df_ab[col+'_b'][missing_vals]
missing_vals = pd.isnull(df[col])
df[col][missing_vals] = df_ab[col+'_a'][missing_vals]
need_to_drop.append('rain_6')
need_to_drop.append('rain_24')
self.df = df.drop(need_to_drop,axis=1)
def _compute_var(self):
geometry = [Point(xy) for xy in zip(self.df.longitude,self.df.latitude)]
crs = {'init':'epsg:4326'}
data = gpd.GeoDataFrame(self.df,crs=crs,geometry=geometry)
GUANGDONG.crs = data.crs
try:
data_guangdong = gpd.sjoin(data,GUANGDONG,how='inner')
self.mean = data_guangdong.mean()
except:
pass
先向大家展示一下读取进来如果不做补缺处理，数据帧是多么凌乱~~~~
image.png
可以说千疮百孔，完全没有办法使用
下面大致介绍这个模块的功能，其实它与序章中提到的CTLReader有类似的结构，先定义一个大类型，然后在每次初始化这种类的时候干以下三件事：
image.png
-
_find_filegroup因为逐小时的观察资料实际上是相当参差不齐的，每小时包含的站点数不同，并且有一些站点当前因各种情况缺值但相邻的前后几个时次都是正常的，考虑到气象数据有一定的连续性，我同时检索了当前以及before和after的一小时的三个文件预备在出现缺值的时候做比较简单的“填补”以确保数据的连续。 -
_read_data读取的前中后三个文件数据后，使用python强大的公式参数化模块lambda x : x[:]+'_a'对相应的数据进行_a或_b的标记，然后用df = pd.concat([data,data_b,data_a],axis=1)对三组数据进行合并,形成如图结构的pandas.DataFrame
| 图1 数据帧的列信息 | 图2 样本列 |
|---|---|
image.png
|
image.png
|
可以看出融合后的数据集存在很多缺测值这主要是因为不同时刻观测站点不同引起的。经过观察可以看出，在当前时刻（没有下标的列），如果有这个站点那么他的经纬度信息始终是完整的，所以我使用df=df[:][pd.notnull(df['longitude'])]提取了当前时刻具有完整经度的站点子集，之后做的工作就是假设同一个站点，如果当前时刻数据缺失，但是它的前一个时刻或后一个时刻具有有效观测则将数据的空缺补上。
image.png
这里优先使用了before的数据，再仍然缺失的情况下则考虑使用after的数据。在补完操作后扔掉多余的数据self.df = df.drop(need_to_drop,axis=1)。
搞定后的数据集结构如下：
| 图1 补缺后的数据帧信息 | 图2 帧实例 |
|---|---|
image.png
|
image.png
|
值得注意的是，虽然我们补缺了很多当前时刻的缺值，但是从前后两次数据帧的信息来看，并没有使得当前总的有效站点数扩大，说明补缺处理基本是成功的（最终的有效站点维持在80个左右，并没有因为引入了前后时次信息而扩大）。
-
_compute_var前一步得到的数据集虽然基本算是丰富且较少缺值了，但是并不是最终我们需要的区域平均。这一步我使用了geopandas模块的交集抓取方案，将点阵的经纬度信息直接投身到GIS系统生成的地理信息shapefile文件上，然后data_guangdong = gpd.sjoin(data,GUANGDONG,how='inner')取了两者的交集，于是新的点阵图像如下：
| 图1 数据详情 | 图2 对应地图的位置 |
|---|---|
image.png
|
image.png
|
可以看出最后交集在广东省内的有意义的测站只有36个，虽然不多但是分布正确而且这只是实例时次，实际观测总是有高有低的。
然后不用管是否有具体意义，交给python直接用data_guangdong.mean()求区域平均，里面就包含了我需要的小时气象要素广东区域平均的观察数据。
image.png
红框内信息代表一个样本的区域平均观测结果，本章之后主要的观测样本就是这些东西。
二. 对模式产品读取的说明
这一部分主要用到了序章中提到的CTLReader但是因为需要根据广东的经纬度信息进行提取，需要修改。
首先考虑到模式产品的经纬度信息一般不会发生变化，并且相当庞大，本例中单单一层的总格点数就在40万左右，这还每算上垂直的分层，所以每一次获取预报值的时候都处理经纬度信息是相当不智的。于是我设计了如下的主程序代码段，用来一次性提取模式产品经纬度信息。
def forcast_in_guangdong(ctl,file):
fcst = CTLReader.CTLReader(ctl,file)
lon,lat=np.meshgrid(fcst.variables['longitude'][:].flatten(),
fcst.variables['latitude'][:].flatten())
res = {}
res['latitude'] = lat.flatten()
res['longitude']= lon.flatten()
res=pd.DataFrame(res)
geometry = [Point(xy) for xy in zip(res.longitude,res.latitude)]
crs = {'init':'epsg:4326'}
data = gpd.GeoDataFrame(res,crs=crs,geometry=geometry)
OBSReader.GUANGDONG.crs = data.crs
datainner=gpd.sjoin(data,OBSReader.GUANGDONG,how='inner',op='within')
return datainner.drop(['index_right','2000','DZZSTLSMJ'
,'QLDJ','QLXS','XZQMC'],axis=1)
这段程序的返回值是去掉了无用信息且保持了格点唯一编号的广东省内经纬度信息。
| 图1 数据详情 | 图2 对应点于地图上 | 图2 放大后珠江口位子 |
|---|---|---|
image.png
|
image.png
|
image.png
|
可以看出虽然只有3km的分辨率，模式的格点数量也远远高于观测，大出去好几个数量级。因而我们取得的这个模式产品的区域平均相比观测的区域平均更具有代表性，两者实际上相比的意义有限，这是一个亟待解决的大命题，这里就不展开讲了，为了专注于本章内容，这里权且当作他们俩有一定关系~~~~
由1W7个点组成。这里的行编号是唯每个点唯一的ID，因此可以利用这个index与数据点阵取（40W点左右）交集方便的得到我需要的广东省内的数据点(这就是pandas强大的检索功能了，40W个点的检索几乎瞬间完成)。
由主程序中这么一个套在CTLReader外面的方法来实现上面的描述，交集取关注信息：
# lonlat_info 为forcast_in_guangdong执行后的带的经纬度信息对象
def forcast_to_result(ctl,file,lonlat_info):
fcst = CTLReader.CTLReader(ctl,file)
res = {}
for var in need_to_read:
res[var]=fcst.variables[var][:].flatten()
res=pd.concat((lonlat_info,pd.DataFrame(res)),axis=1).dropna()
return res.mean()
三. 主程序执行的流程
具体的实现流程包括如下几个部分：
1.从观测站点文件提取观测信息（主要模块为OBSReader）。
2.从模式预报的产品文件中提取主要关注的预报量（主要模块为序章中提到的CTLReader）。
3.数据处理以及为资料做统一的资料清理工作（往往是机器学习模型成功建立最关键的一步）。
4.对清理完成的资料集做分train,test两组。
5.将train组数据“喂”给模型，完成模型训练。
6.用所得模型对test组的资料做“模拟预报”，以验证模型的正确性。
并且由于选用的机器学习方法比较基础，所以上面的第3步尤其复杂，请各位看官重点留意～
1.从观测站点文件提取观测信息
我们首先将每组观测或预报结果看成一个记录（record），为了让记录具有唯一清晰的识别ID，即告诉python如果我给出了一个特定的ID，那么我指的是它对应时刻的观测信息以及预报信息组成的唯一对应的数据序列。根据这个特征我选择了时间作为数据集的唯一ID。
#获取观测结果
features = ['date','temperatura','dewp_temperatura',
'pressure','psl','humidity','u_wind','v_wind']
ObsSummary = namedtuple("ObsSummary",features)
命名了观测数据集的具体元组名字（ObsSummary），它们包含了时间，温度，露点温度，地面气压，海平面气压，U风速，V风速。
path = os.path.abspath('D:\data_synop04')
filelist = os.listdir(path)
filelist.sort()
records = extract_obs_data(path,filelist)
获取并排序观测数据文件，并执行extract_obs_data获取所以列表内文件的记录，统一保存于records。extract_obs_data方法如下：
def extract_obs_data(path,filelist):
records = []
for file in filelist:
obs = OBSReader.obs_from_file(path+os.sep+file)
#这里直接使用OBSReader里定义好的正则搜索对象来检索文件名获取有用的信息
otime = OBSReader.P.search(file).group(2)
time = dateparser.parse("20%s/%s/%s-%s:00:00"%(otime[:2],otime[2:4],otime[4:6],otime[6:8]))
if hasattr(obs,'mean'):
records.append(ObsSummary(
date=time,
temperatura=obs.mean['temperatura'],
dewp_temperatura=obs.mean['dewp_temperatura'],
pressure=obs.mean['pressure'],
psl=obs.mean['psl'],
humidity=obs.mean['humidity'],
u_wind=obs.mean['u_wind'],
v_wind=obs.mean['v_wind']))
return records
调用OBSReader并利用正则表达搜索出的时间字符串调用dateparser.parse生成时间类别对象，然后将obs对象中的区域平均信息直接传递给元组ObsSummary形成一个记录。
重复进行循环，得到的记录组是这样的元组列表(records):
image.png
因为观测资料是文本格式的，而且处理程序写得相对复杂，所以这一步会执行比较久的时间，为了方便重复实验避免不必要的重复读取数据，我这里引入了一个python中强大的中间格式存储库pickle，这是”韩国泡菜“的英文单词，如单词的字面意思，这个模块能够将几乎所有的python对象重内存中”挖“出来，然后"塞入泡菜坛子"（一个较小的二级制专业中间文件）。之后再需要使用的时候只有遵循统一的简单操作就能得到和第一次运行程序得到的对象一摸一样的东西。
请注意，我这里用的是”对象“这个词，代表着这个模块实现的功能不只是将数组写入文件这样简单的操作，甚至还可以存储训练好的模型，有实例定义的方法等等，真的是一个泛用性非常强的中间格式存储库，但如果你追求绝对的安全还是把重要数据写成数据文件保存更好>,<
马上来对刚才生成的records进行泡菜存储吧 ^ . ^
import pickle
with open('records.pkl','wb') as fp:
pickle.dump(records,fp)
#在执行目录会生成一个叫”records.pkl“的文件，完全不用在意里面是什么格式，什么东西，大约就是一缸泡菜
之后想要取出泡菜里的东西只需要执行下面这段简单的操作，里面的东西就能变得和你执行一大段代码获得的全新对象一摸一样。
with open('records.pkl','rb') as fp:
records=pickle.load(fp)
获取到的记录列表对象可以用pandas模块方便的转化成漂亮的数据帧，方便人理解，更方便机器处理。
df = pd.DataFrame(records,columns=features).set_index('date')
得到的pandas.DataFrame实例对象如下：
| 图1 数据帧简单信息 | 图2 数据帧具体内容 |
|---|---|
image.png
|
image.png
|
可以看出，数据的缺值情况还是有少量存在。之后可以通过数据清理解决。
2.从模式预报的产品文件中提取主要关注的预报量
接下来通过预报信息命名为类似的命名元组（FcstSummary）
fcst_features = ['date','CombinedRadarReflective','SurfacePressure',
'SeaLevelPressure','RelativeHumidity2m','Temperature2m',
'SurfaceTemperature','U_wind10m','V_wind10m',
'CombinedRadarReflective_prefcst','SurfacePressure_prefcst',
'SeaLevelPressure_prefcst','RelativeHumidity2m_prefcst','Temperature2m_prefcst',
'SurfaceTemperature_prefcst','U_wind10m_prefcst','V_wind10m_prefcst']
FcstSummary = namedtuple("FcstSummary",fcst_features)
然后以一定方法获取想要的预报信息记录组，为了使得程序可读性更好，这里说明的循环模块把所有的获取动作统一替换为print的输出结果方便显示流程，加强理解，原始代码可以直接参考主程序。
fcst_records=[]
for date in df.index:
validtime = date.to_pydatetime()
fcst_date_before=date.date() - forcast_d
fcst_date=dateparser.parse(date.date().strftime("%Y/%m/%d-%H:%M:%S"))
fcst_date_before=dateparser.parse(fcst_date_before.strftime("%Y/%m/%d-%H:%M:%S"))
timedelta = validtime - fcst_date
timedelta_before = validtime - fcst_date_before
timedelta=int(timedelta.total_seconds()/3600)
timedelta_before=int(timedelta_before.total_seconds()/3600)
filename1=fcstpath+os.sep+"%s%03d"%(fcst_date.strftime("%Y%m%d%H"+os.sep+"mars3km%y%m%d%H"),timedelta)
filename2=fcstpath+os.sep+"%s%03d"%(fcst_date_before.strftime("%Y%m%d%H"+os.sep+"mars3km%y%m%d%H"),timedelta_before)
print('-----------------获取预报信息-----------------')
if os.path.isfile(filename1) and os.path.isfile(filename2):
print('判断文件是否存在')
print('-----------------验证时间为：',validtime,'-----------------')
print('这里调用forcast_to_result方法从',filename1,'获取预报信息，表示离验证时间00-24小时内的预报结果')
#f1=forcast_to_result(ctl,filename1,lonlat_in_guangdong)
print('这里调用forcast_to_result方法从',filename2,'获取预报信息，表示离验证时间24-48小时内的预报结果')
#f2=forcast_to_result(ctl,filename2,lonlat_in_guangdong)
print('将获取的两次预报信息写成以时间为ID的唯一记录')
"""
fcst_records.append(FcstSummary(
date =validtime,
CombinedRadarReflective=f1.cr,
SurfacePressure =f1.ps,
SeaLevelPressure =f1.psl,
RelativeHumidity2m =f1.rh2m,
Temperature2m =f1.t2m,
SurfaceTemperature =f1.ts,
U_wind10m =f1.u10m,
V_wind10m =f1.v10m,
CombinedRadarReflective_prefcst=f2.cr,
SurfacePressure_prefcst =f2.ps,
SeaLevelPressure_prefcst =f2.psl,
RelativeHumidity2m_prefcst =f2.rh2m,
Temperature2m_prefcst =f2.t2m,
SurfaceTemperature_prefcst =f2.ts,
U_wind10m_prefcst =f2.u10m,
V_wind10m_prefcst =f2.v10m,
))
"""
print('-----------------一次操作结束-----------------')
模拟执行以后的输出流程如下：
image.png
通过几个标注时间的检查，可以看出程序稳定的找到了每一个观测时间点对应的预报场文件，并且因为书写时用的是统一系统操作模块import os，代码是可以在任何系统不做任何修改的流畅兼容的。
接下来还是制作预报泡菜的操作，不在赘述：
image.png
3.数据处理以及为资料做统一的资料清理工作
接下来将预报信息记录组转成数据帧，并将名字提取成列表备用，之后合并观测与预报两个数据帧，形成新的总数据帧：
| 代码段 | 总数据帧描述信息 | 数据帧内容展示(部分) |
|---|---|---|
image.png
|
image.png
|
image.png
|
接下来需要介绍下面的概念，首先我们这次关注的温度预报，其实为整个气象业务里相当理想的情况了，因为温度这个量具有很多"friendly"的统计特征，包括是连续变化的要素，变化相对缓慢，统计分布基本服从正态分布等等，这样就使得本章面对具体任务大大理解化了，所以温度这个量无论从什么角度出发都是气象业务里相对容易预报的，换句话说是”用户友好型物理量“。
利用早就被前人发现的温度连续变化的物理特征，我用了下面的方法扩展了每个记录实际预报的信息，既将每个记录点的向前个时次的预报信息归并到一起，这样一个记录就同时拥有了 ： 观测，前一小时的观测，前两小时的观测，前三小时的观测，预报，前一小时的预报，前两小时的预报，前三小时的预报 。这样的9大类信息，之后再慢慢分析这些信息源的可预报关系。
实现代码如下：
def derive_ntime_feature(df,feature,N):
rows = df.shape[0]
ntime_prior_measurements = [None]*N + [df[feature][i-N] for i in range(N,rows)]
col_name = "{}_{}".format(feature,N)
df[col_name] = ntime_prior_measurements
for feature in features_all:
if feature != 'date':
for N in range(1,4):
derive_ntime_feature(df_all,feature,N)
执行后带有标记的数据帧展示如下，边界上的几个时次直接用了缺值替换掉：
| 总数据帧描述信息 | 数据帧内容展示(部分) |
|---|---|
image.png
|
image.png
|
接下来又是显示python强大的列表操作方法的时候了，因为本章只把温度当成预报对象，所以观测里不带有下划线的变量都可以去掉，以下代表段表示这么一个意思：列出一个名叫to_remove的列表，它的内容是一些feature，这些feature应该包含在features_obs里面并且不包含变量temperatura。
to_remove = [feature for feature in features_obs if feature not in ['temperatura']]`
还有另一个剔除的原因是如果你知道了当前时刻的这些观测量，那么就代表你也同时知道了当前时刻的温度，这样就失去了训练模型的意义了。我就得到了一个需要移除的列表名字，同时保留了关注的部分。
image.png
然后利用类似的列表操作方式将df_all里面不需要的部分排除得到：
image.png
利用.describe()方法可以粗略的看出df_all的统计信息，利用简单的分析可以大致看到变量信息是否分布上有严重离群的情况出现：
image.png
可以看出大部分预报及观测的分布没有严重离群，但是纬向分的观测与模式的相对湿度除外，这其实还比较合理，因为纬向风有盛行风向，导致分布情况不均匀有偏态。如图：
image.png
分布上说明4月广东地区平均分纬向以弱的偏东风为主，这很符合实际。
同时可以看到模式的相对湿度预报有一个最大阈值99.144684，同时存在一些奇异点，相对湿度为0，而且有一个有趣的现象，既模式24小时以内的预报均有这种奇异点发生，而24-48小时的预报结果就完全没有这样的事情，示例如图：
image.png
这不符合物理现实，其实是因为模式的每一个起报时次的预报产品水物质含量为0（相对湿度为0），这又牵涉到另一个超大课题模式的起转问题了。。。。，但这种情况处理起来比较麻烦，所以为了方便描述，之后的处理中打算直接扔掉相对湿度和风速变量。
插入一段对以下计算的相关性的统计学介绍：
首先相关系数表示两个量直接的相关性，值域在-1到+1之间，负值表示两个量的变化是反向的既反相关，正值表示两个量的变化趋势一致，绝对值越靠近1则表示两个量之间的相关性越强，绝对值越靠近0则代表两个变量直接的相关性越弱。但是统计学上并没有严格区分哪个相关性范围的值适合作为预报量，按照下表我自行定义了一个大致的区域：
| 相关系数绝对值 | 代表意义 |
|---|---|
| 0.8~1.0 | 非常强的相关 |
| 0.6~0.8 | 强的相关性 |
| 0.4~0.6 | 中等相关性 |
| 0.2~0.4 | 弱相关性 |
| 0.0~0.2 | 非常弱的相关性 |
接下来就该考察各个变量与我关心的temperature之间的关系了(这里发现了一个错别字。。。我把温度写成了temperatura........>.<，之后的程序统一改回来了，之前的懒得改了。。。)
df_all = df_all.dropna()
corr=df_all.corr()[['temperature']].sort_values('temperature')
corr=corr[abs(corr)>0.6]
corr=corr.dropna()
利用以上代码我们可以得到df_all里各种变量对应于温度的相关系数绝对值大于0.6的统计结果，如图所示：
image.png
主要的模型训练量就从这里面挑选了，虽然相关系数绝对值在0.6以下的变量是物理上具有意义的预报量，但这里我们专注于训练线性回归模型，为了方便模型的训练舍弃掉那一部分。仔细比较相关系数的值可以看出，temperature_1/2/3相对于temperature都有很高的相关性，分别是0.97，0.91，0.82，这是符合物理常识的，因为温度是逐渐变化的，如果你知道了前三小时的三次观测温度，那么当前时次的温度几乎可以肯定就在这几个值附近，但这样做对模式训练的意义就大打折扣了，因为如果预报模型需要提供给最近三小时的观测，那么这个时效就太短了不符合需求。除此之外，其他变量有强相关的基本上全是模式预报的变量，因此可以考虑完全使用模式预报的变量训练预报模型。
继续排除掉观测量本身及前几个时次的观测后，得到了大致的模型预报变量的名称列表：
predictors = list(corr.drop(['temperature','v_wind_1','temperature_3','temperature_2','temperature_1'],axis=0).T)
接下来要进行的就是筛选独立变量的工作了，因为线性回归中有一条基本的假设：一个理想的状态是回归模型的各个自变量相互独立，其各自的变量只会导致因变量的变化，因此自变量之间如果存中相关性就会破坏回归模型的稳定性。但实际上这种完全独立的情况是很难达到的，一般的做法是给定一个很小的概率阈值（一般是0.05）进行统计检验，如果自变量之间的相关超过了自变量与因变量的相关的发生概率大于这个阈值则认为这个变量不独立，做剔除处理。这里又会继续涉及到统计学的第一类错误与第二类错误。。。。有兴趣的自己谷歌吧。
交给python去做统计检验其实很简单，利用statsmodels.api这个标准的统计学库就可以方便的进行(这里可能会出现警告，但不影响检验，忽略之)。执行以下操作将 df2 分解成自变量部分与因变量部分，并对自变量部分添加一个常数列，这样就可以用标准化的检验方法开始检验了。
df2 = df_all[['temperature'] + predictors]
df2 = df2.dropna()
import statsmodels.api as sm
X = df2[predictors]
y = df2['temperature']
X = sm.add_constant(X)
X.iloc[:5,:5]
前几个行列的展示如下：
image.png
利用以下代码可以方便的得到当前的统计检验结果报告，如下图：
model = sm.OLS(y,X).fit()
model.summary()
image.png
看起来复杂，但本章的需求只需要重点关注框选部分：
蓝框部分的R-squared（又可以称作coefficient of determination）既决定系数，表征依变数Y的变异中有多少百分比可由控制的自变数X来解释（这里表示解释了93.2%的变化）。另一个信息是修正后的决定系数（Adjusted R-Squared），这样里需要一点点统计学知识，我们知道在其他变量不变的情况下，引入新的变量，总能提高模型的R2。修正R2就是相当于给变量的个数加惩罚项，如果两个模型，样本数一样，R2一样，那么从修正R2的角度看，使用变量个数少的那个模型更优。换句话说，如果两个模型，样本数一样，R2一样，那么从修正R2的角度看，使用变量个数少的那个模型更优。
红框部分，就是上面我提到的阈值判断，如果一个自变量的假设检验概率大于0.05则表示这个模型中这个变量与其他变量是相关的应该剔除，但需要注意每次剔除一个变量则应该看作是一个全新的模型，这样每次只剔除最大的超过0.05的变量，再次做统计检验直到没有变量的假设检验概率大于0.05。可以编写以下的循环来实现：
vf = True
while vf is True:
model = sm.OLS(y,X).fit()
pmax = model.pvalues.max()
if pmax > 0.05:
X=X.drop(model.pvalues.idxmax(),axis=1)
else:
vf = False
最后得到了13个独立性较好的预报变量，终于可以将这些变量当作模型的自变量进行训练了。
image.png
值得注意的是，蓝框中的统计量之差减少了，这表明刚才的筛选不仅减少了模型的自变量数目，也使得模型出现过度拟合（overfitting）的几率减小了。
4.对清理完成的资料集做分A,B两
首先说明为了避免“幸存者偏差”使得有限的训练样本尽量代表广泛一些的信息，在进行分组时需要进行随机抽取。因此我使用了标准的随机抽取方法:
from sklearn.model_selection import train_test_split
X = X.drop('const',axis=1)
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=12)
这里的random_state属性只要给定一个常数，就会让每一次执行train_test_split方法的结果保持一致，方便进行重复测试。test_size=0.2则表示测试用的样本占总样本的20%。
5.将train组数据“喂”给模型，完成模型训练
利用sklearn库里面标准的线性模型进行训练：
from sklearn.linear_model import LinearRegression
regressor = LinearRegression()
regressor.fit(X_train,y_train)
执行完成后就得到了一个名叫regressor的模型对象，这个就是最后的训练结果。你甚至还可以使用pickle库将这个线性模型当作泡菜保存起来，下次接着训练或是使用。
6.用所得模型对test组的资料做“模拟预报”，以验证模型的正确性
预报的实际执行只需要下面这一小步，即可得到预报结果。
prediction = regressor.predict(X_test)
为了方便的分析预报结果，我下面将通过训练模型得到的预报结果称为"end_var"然后将前面通过分组得到的每个时次的两次单纯模式预报结果称为“'Temperature2m','Temperature2m_prefcst”，那么单纯的当前时刻温度观测的结果既为"temperature"。
res=pd.concat([y_test,pd.DataFrame({'end_var':prediction},index=y_test.index),only_fcst],axis=1).dropna()
拼接后的分析数据大约长这样~：
image.png
将上面的数据帧画成折线图可以看得更明显。
image.png
可以看出在测试数据的绝大部分温度转折区域，训练后的模型预报（橙色）明显要优于两次单纯的模式预报结果（绿色与红色），途中标记的箭头大体能看出线性回归方法预报的好（红）坏（黑），通过下面这一段统计信息可以对比总体的差异
from sklearn.metrics import mean_absolute_error,median_absolute_error
print('线性方案的平均绝对误差：',mean_absolute_error(res['temperature'],res['end var']),
'线性方案的绝对误差中位数：',median_absolute_error(res['temperature'],res['end var']))
print('纯3km数值模式的24小时平均绝对误差：',mean_absolute_error(res['temperature'],res['Temperature2m']),
'纯3km数值模式的24小时绝对误差中位数：',median_absolute_error(res['temperature'],res['Temperature2m']))
print('纯3km数值模式的48小时平均绝对误差：',mean_absolute_error(res['temperature'],res['Temperature2m_prefcst']),
'纯3km数值模式的48小时绝对误差中位数：',median_absolute_error(res['temperature'],res['Temperature2m_prefcst']))
结果为：
image.png
最后按照前面得到的回归模型，不再进行重新训练，额外引入了2017年5月份的模式预报结果进行验证得到的结果如下：
image.png
可以看出红色与蓝色线段明显靠近了。
但是效果比test试验要略差，说明随着季节气候的改变等因素，统计方法得到的订正模型正在逐渐退化，因此需要继续发展可以实时更新的预报模型使得模型的效果保持相对稳定。
同时，通过对比前后两次测试的总体结果数据：
image.png
可以看出，数值模式的统计误差是基本固定的。因此是完全有可能通过引入自动更新的机器学习方法达到稳定的降低一定模式系统误差的。
四. 结论与讨论
在本章中我介绍了如何用简单的线性回归算法去预报气温，虽然本章仅仅讨论了一个高度简化环境下，最原始的机器学习算法预报效果，但能看出以下几条：
- 模式预报并不完美，其效果随着时效延长而逐渐变差。
- 机器学习及人工智能方法可以通过大量数据训练纠正一部分模式预报的系统性误差。
- 虽然对数据的清理工作看似繁琐，但核心思想只有一条，对复杂的多种类数据进行标记（addressing），本例子中用的是时间ID，让其拥有唯一确定的ID，以便建立关系数据帧。
- 这仅仅是一个开头，之后有时间再继续介绍其他算法的应用。
- 下一步的打算是使用其他算法，做出空间分布的预报模型而非本章这种单纯的空间平均结果，并试试看一些更有挑战的物理量，比如降水，对流高度，雷达反射率等等。
最后的最后顺便帮同学打个广告，我码字这么轻松就能写这么多主要是多亏了有“航天枸杞”保驾护航~~~_
image.png
image.png
网友评论