置信度检验
<p>[TOC]</p>
<pre><code class="language-python">%matplotlib inline
%load_ext autoreload
%autoreload 2
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>
<pre><code class="language-python">sta_all = pd.read_hdf(r&quot;H:\test_data\input\mpd\Example_data\sta_all.h5&quot;)
print(sta_all)</code></pre>
<pre><code> level time dtime id lon lat OBS ECMWF \
0 0.0 2019-07-01 08:00:00 0 54398 116.6 40.1 25.8 24.724
1 0.0 2019-07-01 08:00:00 0 54410 116.1 40.6 18.9 20.284
5 0.0 2019-07-01 08:00:00 0 54412 116.6 40.7 19.7 19.820
2 0.0 2019-07-01 08:00:00 0 54416 116.9 40.4 25.1 22.864
3 0.0 2019-07-01 08:00:00 0 54419 116.6 40.4 27.5 21.796
... ... ... ... ... ... ... ... ...
9187 0.0 2019-07-31 20:00:00 72 54410 116.1 40.6 18.8 21.556
9191 0.0 2019-07-31 20:00:00 72 54412 116.6 40.7 25.4 21.640
9188 0.0 2019-07-31 20:00:00 72 54416 116.9 40.4 26.7 24.960
9189 0.0 2019-07-31 20:00:00 72 54419 116.6 40.4 27.0 22.920
9190 0.0 2019-07-31 20:00:00 72 54499 116.2 40.2 29.3 23.056
GRAPES
0 23.844
1 21.740
5 21.900
2 22.704
3 23.136
... ...
9187 19.968
9191 21.232
9188 23.404
9189 23.056
9190 22.568
[9192 rows x 9 columns]</code></pre>
<h1>相对技巧置信度</h1>
<p><font face="黑体" color=blue size = 3>score_compare(sta_ob_and_fos0,method,grade_list = None,compare = ">=",s = None,
B = 10000 ,c = 0.95
,ncol = 1 ,title_list = None ,vmax = None ,vmin = None ,vmax_delta = None,
sparsify_xticks = 1, sup_fontsize = 12 ,xlabel =None ,ylabel = "Score" ,
ylabel_delta = "Score Delta" ,width = 10,height = 6,
dpi=300 ,legend_loc = "upper center" ,wspace = None,
hspace = None,save_path = None,show =False)</font> </p>
<p>采用如下方法绘制相对技巧置信图： </p>
<ol>
<li>根据输入的站点数据和检验方法，选择部分数据，将选择的数据进行按照时效进行分组，计算的不同时效的评分。</li>
<li>选取1个时效的数据，计算不同起报时刻（例如365天）的评分，得到2×365个评分</li>
<li>用第2个模式的评分减第1个模式的平方，得到365个相当技巧</li>
<li>采用又放回重采样方法从365个相当技巧值选取365个值（部分重复），取平均值</li>
<li>重复B次(例如10000）步骤4，得到B个平均相当技巧</li>
<li>将步骤5的结果从大到小排序，选取排序在前a%和后a%的位置的两个取值lower 和higher， c = 1-2*a</li>
<li>取(higher - lower) 作为置信区间宽度。</li>
<li>遍历不同时效，重复步骤2-7</li>
<li>绘图检验中间量）</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>sta_ob_a0nd_fos </font></strong></td>
<td style="text-align: left;">实况和预报合并对齐后的数据，形式为站点数据格式如上述例子中的sta_all</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>method </font></strong></td>
<td style="text-align: left;">method中的各类数值型检验指标的函数名称，例如在本页面中已经import meteva.method as mem ，则ts评分的函数名称可以写为mem.ts ，均方根误差计算函数的名称可以写为 mem.rmse。这些函数选项包括<a href="https://www.showdoc.cc/nmc?page_id=2858658548509727">连续型预报</a>，<a href="https://www.showdoc.cc/nmc?page_id=2860336958932349">二分类预报</a>,<a href="https://www.showdoc.cc/nmc?page_id=2859693269266585">多分类预报</a>,<a href="https://www.showdoc.cc/nmc?page_id=3651805946039771">概率预报</a>,<a href="https://www.showdoc.cc/nmc?page_id=3629735872716279">集合预报</a>等类别中的评分</td>
</tr>
<tr>
<td style="text-align: left;"><strong>grade_list</strong></td>
<td style="text-align: left;">该参数用于对连续变量做多种等级阈值的二分类检验，其中包含多个事件是否记录为发生的判断阈值，记其中一个阈值为g，则判断为事件发生的条件是要素值 >= g。该参数缺省时列表中只包含一个取值为1e-30的阈值，由于气象要素精度通常比该缺省值大，因此它相当于将 >0 作为事件发生的判据</td>
</tr>
<tr>
<td style="text-align: left;"><strong>compare</strong></td>
<td style="text-align: left;">比较方法，可选项包括">=",">","<=","<", 是要素原始值和阈值对比的方法，默认方法为">="，即原始值大于等于阈值记为1，否则记为0，compare 为"<="时，原始值小于等于阈值的站点记为1， 这在基于能见度标记大雾事件、低温事件或降温事件等场景中可以用到</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>B</strong></td>
<td style="text-align: left;">bootstrap算法重采样的次数线</td>
</tr>
<tr>
<td style="text-align: left;"><strong>c</strong></td>
<td style="text-align: left;">置信度阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong>ncol</strong></td>
<td style="text-align: left;">当包含多个等级时，子图的列数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>title_list</strong></td>
<td style="text-align: left;">每个子图的标题</td>
</tr>
<tr>
<td style="text-align: left;"><strong>vmax</strong></td>
<td style="text-align: left;">绘制的检验指标图片中y轴坐标的范围的最大值，如果只有一个子图，vmax 应该为一个实数；当有多个子图时，如果vmax为实数，则每个子图y轴的最大值相同，如果vmax为一个实数列表，则每个子图中y轴最大值根据列表中对应值进行设置</td>
</tr>
<tr>
<td style="text-align: left;"><strong>vmin</strong></td>
<td style="text-align: left;">绘制的检验指标图片中y轴坐标的范围的最小值，如果只有一个子图，vmin 应该为一个实数；当有多个子图时，如果vmin为实数，则每个子图y轴的最小值相同，如果vmin为一个实数列表，则每个子图中y轴最小值根据列表中对应值进</td>
</tr>
<tr>
<td style="text-align: left;"><strong>vmax_delta</strong></td>
<td style="text-align: left;">相对技巧的取值范围为 -vmax_delta 至vmax_delta</td>
</tr>
<tr>
<td style="text-align: left;"><strong>spasify_xticks</strong></td>
<td style="text-align: left;">x轴坐标刻度的稀疏倍数，缺省时会不进行稀疏化</td>
</tr>
<tr>
<td style="text-align: left;"><strong>sup_fontsize</strong></td>
<td style="text-align: left;">图片标题的字体大小，其它字体将根据标题字体大小自动设置，其中坐标轴字体大小 = sup_fontsize <em> 0.9, 坐标刻度的字体大小 = sup_fontsize </em> 0.8</td>
</tr>
<tr>
<td style="text-align: left;"><strong>xlabel</strong></td>
<td style="text-align: left;">子图的横坐标名称</td>
</tr>
<tr>
<td style="text-align: left;"><strong>ylabel</strong></td>
<td style="text-align: left;">评分子图的纵坐标名称</td>
</tr>
<tr>
<td style="text-align: left;"><strong>ylabel_delta</strong></td>
<td style="text-align: left;">技巧子图的纵坐标名称</td>
</tr>
<tr>
<td style="text-align: left;"><strong>width</strong></td>
<td style="text-align: left;">图片的宽度，缺省时程序自动设置</td>
</tr>
<tr>
<td style="text-align: left;"><strong>height</strong></td>
<td style="text-align: left;">图片的高度，缺省时程序自动设置</td>
</tr>
<tr>
<td style="text-align: left;"><strong>dpi</strong></td>
<td style="text-align: left;">绘图的dpi参数，用法同matplotlib中dpi参数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>legend_loc</strong></td>
<td style="text-align: left;">legend的位置</td>
</tr>
<tr>
<td style="text-align: left;"><strong>wspace</strong></td>
<td style="text-align: left;">子图之间的横向间隔</td>
</tr>
<tr>
<td style="text-align: left;"><strong>hspace</strong></td>
<td style="text-align: left;">子图之间的纵向间隔</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,则该参数会自动切换为T于1.</td>
</tr>
</tbody>
</table>
<p><strong>调用示例：</strong></p>
<pre><code class="language-python">mpd.score_compare(sta_all,mem.rmse,title_list = [&quot;温度预报RMSE评分和相对大小的置信度检验&quot;],xlabel = &quot;时效(单位:小时)&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=0d57d669c06fe19c07518878d9e312d6&amp;file=file.png" alt="" /></p>
<p>上面的图中下半部分子图中黑边白面柱子的上下跨度就是根据bootstrap获得置信区间大小。</p>
<pre><code class="language-python">mpd.score_compare(sta_all,mem.ts,grade_list = [20,25],title_list = [&quot;温度预报TS评分(&gt;=20℃)&quot;,&quot;温度预报TS评分(&gt;=25℃)&quot;],ncol =2,xlabel = &quot;时效(单位:小时)&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=787888210aa4420dda04ab3e140e03b2&amp;file=file.png" alt="" /></p>
<h1>评分叠加置信区间</h1>
<p><font face="黑体" color=blue size = 3>score_confidence(sta_ob_and_fos0,method,grade_list = None,compare = ">=",s = None,B = 1000 ,c = 0.95,
plot = "bar",**kwargs)</font> </p>
<p>采用如下方法绘制相对技巧置信图： </p>
<ol>
<li>根据输入的站点数据和检验方法，选择部分数据，将选择的数据进行按照时效进行分组，计算的不同时效的评score分。</li>
<li>选1个模式取1个时效的数据，计算不同起报时刻（例如365天）的评分，得×365个评巧</li>
<li>采用又放回重采样方法从365个相当技巧值选取365个值（部分重复），计算得到1个平均值</li>
<li>对步骤3重复B次，得到B个平均值</li>
<li>将步骤4的结果从大到小排序，选取排序在前((1-c)/2)%和后((1-c)/2)%的位置的两个取值lower 和higher</li>
<li>取(higher - lower2 作为置信区间宽。</li>
<li>遍历不同时效，重复步骤2-5</li>
<li>绘图时将置信区间的中心点设置在步骤1计算的到的单时效总体评分score的位置*</li>
</ol>
<p><font face="黑体" color=red size = 3>注意：对于ME\MAE这样的检验指标来说，采用有放回采样取平均值的方法获取的置信区间是合理的。但rmse、ts等指标采用有放回采样取平均值的方法获取的置信区间不太严格，但考虑到非常严格的置信统计较为复杂、计算量也会大增，因此暂时统一采用上面的方案</font> </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_ob_and_fos0 </font></strong></td>
<td style="text-align: left;">实况和预报合并对齐后的数据，形式为站点数据格式如上述例子中的sta_all</td>
</tr>
<tr>
<td style="text-align: left;"><strong><font face="黑体" color=blue size = 5>method </font></strong></td>
<td style="text-align: left;">method中的各类数值型检验指标的函数名称，例如在本页面中已经import meteva.method as mem ，则ts评分的函数名称可以写为mem.ts ，均方根误差计算函数的名称可以写为 mem.rmse。这些函数选项包括<a href="https://www.showdoc.cc/nmc?page_id=2858658548509727">连续型预报</a>，<a href="https://www.showdoc.cc/nmc?page_id=2860336958932349">二分类预报</a>,<a href="https://www.showdoc.cc/nmc?page_id=2859693269266585">多分类预报</a>,<a href="https://www.showdoc.cc/nmc?page_id=3651805946039771">概率预报</a>,<a href="https://www.showdoc.cc/nmc?page_id=3629735872716279">集合预报</a>等类别中的评分函数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>grade_list</strong></td>
<td style="text-align: left;">该参数用于对连续变量做多种等级阈值的二分类检验，其中包含多个事件是否记录为发生的判断阈值，记其中一个阈值为g，则判断为事件发生的条件是要素值 >= g。该参数缺省时列表中只包含一个取值为1e-30的阈值，由于气象要素精度通常比该缺省值大，因此它相当于将 >0 作为事件发生的判据</td>
</tr>
<tr>
<td style="text-align: left;"><strong>compare</strong></td>
<td style="text-align: left;">比较方法，可选项包括">=",">","<=","<", 是要素原始值和阈值对比的方法，默认方法为">="，即原始值大于等于阈值记为1，否则记为0，compare 为"<="时，原始值小于等于阈值的站点记为1， 这在基于能见度标记大雾事件、低温事件或降温事件等场景中可以用到</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>B</strong></td>
<td style="text-align: left;">bootstrap算法重采样的次数</td>
</tr>
<tr>
<td style="text-align: left;"><strong>c</strong></td>
<td style="text-align: left;">置信度阈值</td>
</tr>
<tr>
<td style="text-align: left;"><strong>plot</strong></td>
<td style="text-align: left;">取bar或line，对应bar图和折线图</td>
</tr>
<tr>
<td style="text-align: left;">kwargs</td>
<td style="text-align: left;">检验方法 meteva.method 和meteva.base.bar或meteva.base.plot 中的可选参数，具体用法参见下面的示例</td>
</tr>
</tbody>
</table>
<pre><code class="language-python">mpd.score_confidence(sta_all,mem.rmse,c = 0.8,vmin = 0,plot = &quot;line&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=20fc43ec3ba58b75f8523d40c7b846a6&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">mpd.score_confidence(sta_all,mem.rmse,c = 0.8,vmin = 0,plot = &quot;bar&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=f976cba2cd368e6c734b92ae17e487a4&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">mpd.score_confidence(sta_all,mem.ts,grade_list = [20,25],vmin = 0,plot = &quot;line&quot;)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=3ba3c62b5c1c3c59e0aa6a09f3976032&amp;file=file.png" alt="" /></p>
<pre><code class="language-python">mpd.score_confidence(sta_all,mem.ts,grade_list = [20,25],bar_width = 0.3,vmin = 0)</code></pre>
<p><img src="https://www.showdoc.com.cn/server/api/attachment/visitFile?sign=0ff728454e7b44d5300291a8f8e9f4ac&amp;file=file.png" alt="" /></p>