斯普罗岛最大风速预测
翻译自:http://scipy-lectures.github.com/intro/summary-exercises/stats-interpolate.html
这个练习目的在于预测每50年的最大风速,如果没有测量的话。可以获得的数据仅仅由位于丹麦的斯普罗气象站测量了21年。首先,统计步骤将给出,然后通过从scipy.interpolate模块中的函数演示。最后有兴趣的读者将被邀请用稍微不同的方法从原始数据计算结果。
统计方法
每年最大风速应该符合正态概率密度函数。然而,这个函数将并不会被估计,因为它从最大风速给出概率。这将是分位函数的角色,这个练习的目标将会是找到它。在当前模型,每五十年发生的最大风速被定义为%2上分位。
i
累积概率p_i
被定义为p_i = i/(N+1)
,其中N=21,有测量的年数。因此计算每年最大风速的累积概率将是可能的。从这些实验点,scipy.interpolate模块将在拟合分位函数时非常有用。最后这五十年的最大值将被从上%2分位数的累积概率函数中求得。
计算累积概率
每年的最大风力已经被计算并以numpy格式保存在文件examples/max-speeds.npy中,因此它们将通过numpy载入:
In [1]: import numpy as np
In [2]: max_speeds = np.load('max-speeds.npy')
In [3]: years_nb = max_speeds.shape[0]
p_i
,相应的值将是:
In [5]: cprob = (np.arange(years_nb, dtype=np.float32) + 1) / (years_nb + 1)
它们被假定用来拟合给定风速:
一元样条预测
UnivariateSpline
类来估计,这个类可以代表源于数据点的样条。默认的行为是构建一个3级(degree)的样条,数据点可以根据它们的可靠性有不同的权重。变体有InterpolateUnivariateSpline
和LSQUnivariateSpline
,它们的误差检查将会改变。在需要二维样条的时候,BivariateSpline
类一族被提供。所有这些一维和二维样条使用FITPACK Fortran子程序,这就是为何通过分别代表和估值样条的splrep
和splev
函数,一个低等的库存取可以使用。而且没有使用FITPACK参数的插值函数也被提供给更简单的使用(参见’interp1d’,interp2d
,barycentric_interpolate
等等)。UnivariateSpline
将被使用,因为3级的样条似乎拟合很正确:
In [7]: from scipy.interpolate import UnivariateSpline
In [8]: quantile_func = UnivariateSpline(cprob,sorted_max_speeds)
分位函数想在将从整个概率范围求值:
In [9]: nprob = np.linspace(0, 1, 1e2)
In [10]: fitted_max_speeds = quantile_func(nprob)
%2
在目前的模型中,每五十年发生的最大风速被定义为上%2分位。结果,累积概率值将是:
In [11]: fifty_prob = 1. - 0.02
所以每五十年发生的暴风风速可被猜测,通过:
In [12]: fifty_wind = quantile_func(fifty_prob)
In [13]: fifty_wind
Out[13]: array(32.97989825386221)
结果现在集中在一个Matplotlib图像上:
Gumbell分布练习
感兴趣的读者现在被邀请通过使用二十一年测量的风速来做个练习。这个测量周期是大约90分钟(原始周期大约10分钟,但是文件大小为了让练习更简单已经被减小)。数据以numpy格式存储在文件sprog-windspeeds.npy中。在你完成练习之前,不要看绘图源码。
- 第一步将是使用numpy并绘制matplotlib条形图,找到每年最大值。
- 第二步是使用Gumble分布——其累积概率
p_i
- 被定义为
-log(-log(p_i))
- ,来拟合一个线性分位函数(记住你可以定义
UnivariateSpline
- 的级数(degree))。绘制每年最大风速与Gumble分布应该给出以下图像:
- 最后一步是找到每五十年发生的最大风速34.23m/s