显著性检验
<p>[TOC]</p>
<pre><code class="language-python">import meteva.base as meb
import meteva.method as mem
import meteva.product as mpd
import numpy as np
import datetime
import copy
import matplotlib.pyplot as plt
import pandas as pd</code></pre>
<h1>显著性检验问题引出</h1>
<p>在预报检验中，分析数值模式和客观预报产品的系统性偏差至关重要，因为本质真正的随机误差是无法订正的，无论是预报员的主观订正还是客观算法的订正都是建立在对已有预报产品的系统性偏差的发掘或认识的基础上。 这里所说的系统性偏差并非指一个固定的值，而是可能和季节、区域或天气形势有关系的一个变化量，因此在检验实践中，通常会采用分类检验的策略，例如不同月份的预报分别统计一个误差分布。问题是，分类越小细，每个类别的样本数越少，从中统计出的误差还是显著的吗？</p>
<p>下面是对ECMWF和CMA_GFS模式的2019年7月份的数据进行误差统计的示例，从结果来看模式的预报存在系统性的偏差，但它们显著吗？</p>
<pre><code class="language-python">#读入温度预报的数据
sta_all = pd.read_hdf(r"H:\test_data\input\mpd\Example_data\sta_all.h5")
#对7月份的24小时时效预报样本进行平均误差的统计，绘制误差的空间分布图
result = mpd.score_id(sta_all,mem.me,s = {"dtime":24,"month":7},ncol = 2,sup_fontsize = 8)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=8c6fc2529059f1a1b4b7c7b16b8c5756&file=file.png" alt="" /></p>
<p>仅从上面图形展示的结果来看，ECMWF和GRAPES模式在所分析区域的大部分站点上是负偏差，在区域西北角一个站点上都是正偏差，从偏差幅度来看，负偏差最低达到-4摄氏度，正偏差最高超过2℃。这种幅度能否说明温度预报存在系统性偏差？亦或者它们仅仅是少量样本统计中出现的随机性？毕竟业务中温度预报出现5℃以上偏差都是常有之事。为了回答这个问题，可以使用本模块的功能来进行定量验证。</p>
<h1>区间检测</h1>
<p><font face="黑体" color=Blue size=3><strong>mean_estimate(error_array, numrep=1000,alpha=0.05)</strong></font> </p>
<p>计算误差的平均值，并通过有放回重采样的方法统计误差均值的变化区间，以alpha=0.05为例，具体方法为：</p>
<ol>
<li>通过有放回的重采样方法，从误差序列中采样numrep次</li>
<li>计算N次采用所得误差值的平均，这作为一次随机采样平均值</li>
<li>重复步骤1-2共1000次，得到1000种平均值</li>
<li>将1000个平均值从小到大排序</li>
<li>取1000个序列的2.5%和97.5%分位值作为误差均值的变化区间。</li>
</ol>
<table>
<thead>
<tr>
<th style="text-align: left;">参数</th>
<th style="text-align: left;">说明</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>error_array</font></strong></td>
<td style="text-align: left;">k×n的误差矩阵，n代表站点（或格点）数，k代表每个站点（格点）上的误差时间序列长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong>numrep</strong></td>
<td style="text-align: left;">重采样次数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>alpha</strong></td>
<td style="text-align: left;">显著性水平</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;">DataFrame，包含Lower,Estimate,Upper三列，分别对应每个格点上的numrep次采用得到的均值的2.5%分位值，一般平均值，numrep次采用得到的均值的97.5%分位值</td>
</tr>
</tbody>
</table>
<p>为了说明区间采样的含义，我们以下先通过两个理想示例加以说明</p>
<pre><code class="language-python">t = np.arange(0,6,0.1) #t代表时间序号
obs = np.sin(t) #生成一个代表观测的序列
fst1 = obs + 0.2 #生成具有固定偏差值的一种预报
fst2 = obs + 0.4*np.random.rand(len(t)) #生成具有系统性正偏差叠加一些随机误差的预报
fst3 = obs + 0.1*np.random.randn(len(t)) # 生成个别时刻有较大误差，但其他时刻没有误差的预报
fst3[30] +=12
plt.plot(t,obs,label = "obs") #绘制观测
plt.plot(t,fst1,label = "fst1") # 绘制第一种预报
plt.plot(t,fst2,label = "fst2") #绘制第二种预报
plt.plot(t,fst3,label = "fst3") #绘制第3种预报
plt.legend()</code></pre>
<pre><code><matplotlib.legend.Legend at 0x268bd694448></code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=513f25112906d6ecf348889d3ada0c6d&file=file.png" alt="" /></p>
<pre><code class="language-python">#一下是第一个预报的误差均值和区间
error = fst1 - obs
result= mem.significance.mean_estimate(error)
print(result)</code></pre>
<pre><code> Lower Estimate Upper
0 0.2 0.2 0.2</code></pre>
<p>在上面的示例中，由于误差是固定值0.2，所有按一般方法统计的误差均值和误差均值的变化范围上下限都是0.2</p>
<pre><code class="language-python">#一下是第一个预报的误差均值和区间
error = fst2 - obs
result = mem.significance.mean_estimate(error)
print(result)</code></pre>
<pre><code> Lower Estimate Upper
0 0.202706 0.230356 0.257584</code></pre>
<p>np.random.rand的均值是0.5，因此在上面的示例中于误差是均值为0.2的一个随机序列，按一般方法统计的误差均值大约是0.2。采用重采用方法得到误差均值的变化范围大约是[0.18，0.23]</p>
<pre><code class="language-python">error = fst3 - obs
result = mem.significance.mean_estimate(error)
print(result)</code></pre>
<pre><code> Lower Estimate Upper
0 -0.020486 0.19656 0.602062</code></pre>
<p>np.random.randn的均值是0，因此在上面的示例中于误差序列的大部分是围绕0附近的随机数，但是序列中有一个幅度为12的大误差样本，
按一般方法统计的误差均值是0.2。而随机采样如果没有采中大误差样本，平均值则会在0附近，甚至低于0。</p>
<h1>单点偏差显著性指标</h1>
<p><font face="黑体" color=Blue size=3><strong>local_sig(ob,fo,numrep=1000, alpha=0.05)</strong></font> </p>
<p>基于随机重采用技术统计单点的误差序列的偏差显著性指标
判断单个误差序列存在系统性偏差的方法：
上面三个示例的平均误差都是0.2左右，但它们并非都存在系统性偏差，实际上误差序列1就是有系统性偏差构成的，误差序列12是由0.2的系统性偏差叠加一个随机误差构成，误差序列3是随机误差叠加一个异常值构成，异常值能提升平均值，但是它不是系统性偏差。通过单点偏差显著指标初步判断一个误差序列是否存在系统性偏差：</p>
<ul>
<li>index = |一般平均值| -（区间上界-区间下界）/2 </li>
</ul>
<p>当index小于0时，上面的指标的代表的含义是平均值离原点0的距离小于区间跨度的一半，也就代表误差均值并没有超过随机性的范围，上面的示例3就属于这种情况；反之则表明误差均值超过了随机性的范围，上面的示例1和2属于这种情况。</p>
<p>结合上面的示例重新来审视区间检测为什么要做有放回抽样的平均，目的就是通过随机的方式排除那些频次不高但对均值影响很大的样本，例如示例3中区间下限的计算结果就是排查了个别大幅偏差后的统计值。</p>
<table>
<thead>
<tr>
<th style="text-align: left;">参数</th>
<th style="text-align: left;">说明</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>ob</font></strong></td>
<td style="text-align: left;">k×n的误差矩阵，n代表站点（或格点）数，k代表每个站点（格点）上的误差时间序列长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>fo</font></strong></td>
<td style="text-align: left;">fo比Ob.shape多一维或者保持一致，若多出一维时，fo是m×k×n的误差矩阵，m代表预报成员数，n代表站点（或格点）数，k代表每个站点（格点）上的误差时间序列长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong>numrep</strong></td>
<td style="text-align: left;">重采样次数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>alpha</strong></td>
<td style="text-align: left;">显著性水平</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;">m×n的矩阵，每个元素是一个预报成员的一个站点上偏差显著性指标</td>
</tr>
</tbody>
</table>
<pre><code class="language-python">index = mem.significance.local_sig(obs,fst3)
print(index)</code></pre>
<pre><code>[-0.12005932]</code></pre>
<pre><code class="language-python">fst = np.array([fst1,fst2,fst3])
index = mem.significance.local_sig(obs,fst)
print(index)</code></pre>
<pre><code>[[ 0.2 ]
[ 0.20084421]
[-0.11968074]]</code></pre>
<p><strong>当单点偏差显著性指标小于0，就一定是存在系统性偏差吗？不一定。</strong> 原因包括两个方面：</p>
<ol>
<li>单点偏差显著性指标的过程有随机性，区间上下界就是通过随机又放回抽样方式获得；</li>
<li>可能单点在长时间并不存在系统偏差，但被检验的时间长度有限，样本数目较少，受随机性影响也可能表现为系统性偏大或偏小。</li>
</ol>
<p>对于原因2，我们可以用下面的试验来进行说明。</p>
<pre><code class="language-python"># 生成一组误差序列，时间长度长60，格点（或站点）数为100，每个点上的误差都是随机正态分布的
obs = np.zeros((60,100))
fst = obs + np.random.randn(60,100)
error = fst - obs
#统计每个点上误差序列对应的一般平均值和均值区间
result = mem.significance.mean_estimate(error)
plt.plot(result["Lower"],label="区间下界")
plt.plot(result["Estimate"],label="一般均值")
plt.plot(result["Upper"],label="区间上界")
plt.legend()
plt.hlines(0,100,0)
plt.ylabel("平均误差")
plt.xlabel("格点（站点）编号")</code></pre>
<pre><code>Text(0.5, 0, '格点（站点）编号')</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=cd950b9064cee5e7917f2715a6e412f5&file=file.png" alt="" /></p>
<pre><code class="language-python">index = mem.significance.local_sig(obs,fst)
plt.plot(index)
plt.hlines(0,100,0,color = "k")
plt.ylabel("单点偏差显著指标")
plt.xlabel("格点（站点）编号")</code></pre>
<pre><code>Text(0.5, 0, '格点（站点）编号')</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=fcae86fea7c0be381a21ff3d4aac4b66&file=file.png" alt="" /></p>
<pre><code class="language-python">index_below_0 = index[index>0]
rate = len(index_below_0)/len(index)
print("偏差表现为显著的格点（站点）数比例："+str(rate))</code></pre>
<pre><code>偏差表现为显著的格点（站点）数比例：0.06</code></pre>
<h1>整场显著性</h1>
<p><font face="黑体" color=Blue size=3><strong>field_sig(ob,fo,member_list = None,numrep=1000, alpha=0.05)</strong></font> </p>
<p>判断整场是否为显著</p>
<p>从上面的示例可以看出，即使是随机的误差场，也会有一定比例的格点（站点）上偏差表现为显著。如果误差场真的有系统性偏差，那这个比例应该会更高。现在问题转化为，当偏差表现为显著的格点（站点）数比例小于多少时，可以认为它仅仅是随机造成，达到多少时可以认定它包含了系统性误差？显著性检验算法给出的判别方法如下：</p>
<ol>
<li>用一个长度为k的随机矢量和每个格点（共n个）的时间序列求相关</li>
<li>用t检验方法判断每个相关系数是否达到显著性标准</li>
<li>计算相关系数达到显著的格点数的比例</li>
<li>重复1-3,1000次，得到1000个比例，并从小大排序</li>
<li>取1000个比例95%分位值作为整场显著的比例阈值，当偏差表现为显著的格点（站点）数比例超过阈值时，可以认为误差场存在系统性偏差</li>
</ol>
<table>
<thead>
<tr>
<th style="text-align: left;">参数</th>
<th style="text-align: left;">说明</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>ob</font></strong></td>
<td style="text-align: left;">k×n的误差矩阵，n代表站点（或格点）数，k代表每个站点（格点）上的误差时间序列长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>fo</font></strong></td>
<td style="text-align: left;">fo比Ob.shape多一维或者保持一致，若多出一维时，fo是m×k×n的误差矩阵，m代表预报成员数，n代表站点（或格点）数，k代表每个站点（格点）上的误差时间序列长度</td>
</tr>
<tr>
<td style="text-align: left;"><strong>member_list</strong></td>
<td style="text-align: left;">预报成员的名称列表</td>
</tr>
<tr>
<td style="text-align: left;"><strong>numrep</strong></td>
<td style="text-align: left;">重采样次数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>alpha</strong></td>
<td style="text-align: left;">显著性水平</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;">包含3列数据的DataFrmae，各列分别是偏差表现为显著的格点（站点）数比例和判断整场为显著所需的比例阈值，误差场中是否存在显著性偏差的判断结果</td>
</tr>
</tbody>
</table>
<pre><code class="language-python">obs = np.zeros((60,100))
fst0 = obs + np.random.randn(60,100)+0.0
fst1 = obs + np.random.randn(60,100)+0.1
fst2 = obs + np.random.randn(60,100)+0.2
fst = np.array([fst0,fst1,fst2])
result = mem.significance.field_sig(obs,fst,member_list = ["fst1","fst2","fst3"])
print(result)</code></pre>
<pre><code> sig_rate rate_threshold is_sig
fst1 0.08 0.15 False
fst2 0.16 0.15 True
fst3 0.31 0.15 True</code></pre>
<p>上面的示例中包含三种预报，三种预报包含的系统性偏差幅度逐渐加大，对应的偏差表现为显著的格点（站点）数比例也逐渐由增大，并最终超过了阈值，被算法识别出来。</p>
<p>在整场显著性判别步骤1中为什么用同一个随机矢量同时和所有格点上的时间序列分布求相关？为什么不是用n组不同的随机向量和n个格点上的序列求相关？为什么不是用n组随机向量和另外n组随机向量求相关？
算法之所以这样设计，是因为相邻的格点（站点）场可能存在相似性，可以想象一个极端的情况：100个格点（站点）是完全是相同的随机时间序列，偏差表现为显著的格点（站点）数比例大部分情况下是0，也有一定概率是1。当比例为1时，也并不能判断整场是显著。也就是说当相邻格点之间存在相似性时，判断为整场显著所需的比例需要提高。步骤1的目的就是模拟这种格点相似性带来的影响的方式，当一个格点上相关达到显著，临近区域也会有一大片达到显著，也就是比例很高，这个高比例最终可能影响步骤5中95%分位值，从而抬升比例阈值。
在上面的示例中，每个格点上误差完全是不相关的，对应的比例阈值是0.15左右，下面我们构建一个具有一定相关性的格点场，看一下对比例阈值的影响。</p>
<pre><code class="language-python">obs = np.zeros((60,100))
fst3 = np.zeros((60,100))
fst4 = np.zeros((60,100))
fst5 = np.zeros((60,100))
base_error = np.random.randn(60) # 不同预报的共同扰动部分，该部分占比越大，格点之间相关性越大
for i in range(100):
fst3[:,i] = obs[:,i] + 1 * base_error + np.random.randn(60)
fst4[:,i] = obs[:,i] + 2 * base_error + np.random.randn(60)
fst5[:,i] = obs[:,i] + 2 * base_error + np.random.randn(60) + 0.5
fst = np.array([fst3,fst4,fst5])
result = mem.significance.field_sig(obs,fst,member_list = ["fst3","fst4","fst5"])
print(result)</code></pre>
<pre><code> sig_rate rate_threshold is_sig
fst3 0.01 0.45 False
fst4 0.00 0.66 False
fst5 0.01 0.65 False</code></pre>
<p>从上面的示例看到当格点之间相关性提高的时候，比例阈值会相应提升，但因为它也是基于随机方法生成的，因此比例阈值还是会存在一些随机性。</p>
<h1>显著性检验（基于站点数据）</h1>
<p><font face="黑体" color=Blue size=3><strong>field_sig_id(sta_all,s = None,numrep=1000,alpha=0.05,plot = None,save_path = None,show = False,</strong>kwargs)**</font> </p>
<p>基于pandas格式站点数据开展显著性检验统计。因为显著性检验要求数据是规整的，因此如果某个时刻某些站点上数据缺失了，会以0补齐。</p>
<table>
<thead>
<tr>
<th style="text-align: left;">参数</th>
<th style="text-align: left;">说明</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=Blue size=5>sta_all</font></strong></td>
<td style="text-align: left;">实况和预报合并对齐后的数据，形式为站点数据格式如上述例子中的sta_all</td>
</tr>
<tr>
<td style="text-align: left;"><strong>s</strong></td>
<td style="text-align: left;">用于选择数据样本的字典参数，具体的参数说明可参见meb.sele_by_dict中的<a href="https://www.showdoc.cc/meteva?page_id=3975604785954540"><font face="黑体" color=red size=5>s</font></a>参数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>numrep</strong></td>
<td style="text-align: left;">重采样次数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>alpha</strong></td>
<td style="text-align: left;">显著性水平</td>
</tr>
<tr>
<td style="text-align: left;"><strong>plot</strong></td>
<td style="text-align: left;">该参数为"scatter"时，会将单点偏差显著性指标的空间分布绘制成图片</td>
</tr>
<tr>
<td style="text-align: left;"><strong>save_path</strong></td>
<td style="text-align: left;">该参数不为None时将图片结果输出值save_path</td>
</tr>
<tr>
<td style="text-align: left;"><strong>show</strong></td>
<td style="text-align: left;">该参数不为None时在屏幕显示图片，如果生成了图片save_path又为None,则该参数会自动切换为True</td>
</tr>
<tr>
<td style="text-align: left;"><strong>kwargs</strong></td>
<td style="text-align: left;">绘图函数meb.scatter_sta中的可选参数，具体用法参见下面的示例</td>
</tr>
<tr>
<td style="text-align: left;"><font face="黑体" color=blue size=5>return</font></td>
<td style="text-align: left;">返回一个元组，其包含2个元素，内容分别是整体显著性检验情况和各单点显著性检验指标</td>
</tr>
</tbody>
</table>
<pre><code class="language-python">index_field,index_station = mem.significance.field_sig_id(sta_all,s = {"dtime":24})
print("整场偏差显著性情况：")
print(index_field)
print("-----------------")
print("各单点偏差显著性指标")
print(index_station)</code></pre>
<pre><code>整场偏差显著性情况：
sig_rate rate_threshold is_sig
ECMWF 1.000000 0.666667 True
GRAPES 0.833333 0.666667 True
-----------------
各单点偏差显著性指标
level time dtime id lon lat ECMWF GRAPES
2577 0 2019-07-01 08:00:00 24 54398 116.6 40.1 1.748240 -0.600413
2578 0 2019-07-01 08:00:00 24 54410 116.1 40.6 2.132018 4.568378
2582 0 2019-07-01 08:00:00 24 54412 116.6 40.7 0.632596 0.419011
2579 0 2019-07-01 08:00:00 24 54416 116.9 40.4 4.784124 2.253235
2580 0 2019-07-01 08:00:00 24 54419 116.6 40.4 3.563565 1.201484
2581 0 2019-07-01 08:00:00 24 54499 116.2 40.2 5.155165 2.737828</code></pre>
<p>上面的统计结果表明，ECMWF和GRAPES模式在分析区域内是都是存在系统性误差的，回应本模块一开头提出的问题。</p>
<pre><code class="language-python">#统计显著性，并将各单点偏差显著性指标绘制成图片
result = mem.significance.field_sig_id(sta_all,s = {"dtime":24},plot = "scatter",ncol = 2)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=bd6c96c0c2f973f34776fcbcab85c1b3&file=file.png" alt="" /></p>
<p>本模块现有的功能还没有考虑到误差序列中，可能存在个别持续了一段时间且偏差很大的过程，它可能让整个误差序列显得有系统性偏差，但实际上一次过程也有很大的偶然性，为了回应这类问题Elmore提出了解决方案，具体方法见下面的参考文献。</p>
<p><strong>参考文献</strong>：
[1] Elmore K L , Baldwin M E , Schultz D M . Field Significance Revisited: Spatial Bias Errors in Forecasts as Applied to the Eta Model[J]. Monthly Weather Review, 1913, 134(2):519-531.</p>