斯普罗岛最大风速预测

翻译自: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)的样条,数据点可以根据它们的可靠性有不同的权重。变体有InterpolateUnivariateSplineLSQUnivariateSpline,它们的误差检查将会改变。在需要二维样条的时候,BivariateSpline类一族被提供。所有这些一维和二维样条使用FITPACK Fortran子程序,这就是为何通过分别代表和估值样条的splrepsplev函数,一个低等的库存取可以使用。而且没有使用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图像上:


source code

Gumbell分布练习

感兴趣的读者现在被邀请通过使用二十一年测量的风速来做个练习。这个测量周期是大约90分钟(原始周期大约10分钟,但是文件大小为了让练习更简单已经被减小)。数据以numpy格式存储在文件sprog-windspeeds.npy中。在你完成练习之前,不要看绘图源码。

  • 第一步将是使用numpy并绘制matplotlib条形图,找到每年最大值。


source code

  • 第二步是使用Gumble分布——其累积概率

p_i

  • 被定义为

-log(-log(p_i))

  • ,来拟合一个线性分位函数(记住你可以定义

UnivariateSpline

  • 的级数(degree))。绘制每年最大风速与Gumble分布应该给出以下图像:


source code

  • 最后一步是找到每五十年发生的最大风速34.23m/s