Bootstrap置信区间和GEV拟合pdf

** Bootstrap置信区间和GEV拟合pdf **

1. 置信区间

置信区间是总体参数估计的一个界限，用于量化估计的不确定性。另外，置信区间是一个范围的可能性。 真正的模型性能可能在这个范围之外。

1.1 分类精度的置信区间

如果给定输入数据，预测它们的标签，通常用分类准确率（accuracy）或分类误差（Error，与准确率相反）来描述分类预测模型的性能，分类准确率或分类误差是一个比例，别名：伯努利审判(Bernoulli trial)。eg: 董某人用Wrf模式模拟了10次兰州沙尘过程，但是只有7次模拟成功，那么模型的分类准确率为70%。分类误差区间半径计算公式：

interval = z * sqrt( (error * (1 - error)) / n)

分类准确率区间半径计算公式：

interval = z * sqrt( (accuracy * (1 - accuracy)) / n)

公式中的interval是置信区间的半径，error和accuracy是分类误差和分类准确率，n是样本大小，sqrt是平方根函数，z是高斯分布的临界值。用术语表述，这就是二项式比例置信区间。

1.2 非参数置信区间

如果我们不知道性能指标的分布情况或者不知道计算置信区间的具体方法或者所拥有数据量太少，在这些情况下我们可以采用 bootstrap重采样方法计算置信区间。 任意总体统计的置信区间都可以用***bootstrap***以一种分布无关法（distribution-free）进行估计。 bootstrap是一种模拟蒙特卡罗方法，其中样本是从固定的有限数据集中有放回的抽取出来的，并且在每个样本上估计一个参数。
Python代码实现：

```python
import numpy as np
def average(data):
return sum(data) / len(data)
def bootstrap(data, B, c, func):
**#计算bootstrap置信区间
#:param data: array 保存样本数据
# :param B: 抽样次数 通常B>=1000
#:param c: 置信水平
#:param func: 样本估计量
# :retrn: bootstrap置信区间上下限
array = np.array(data)
#将数据赋值到np的array数组里
n = len(array)
#数据长度
sample_result_arr = []
#建立空数组，动态数组
for i in range(B):
index_arr = np.random.randint(0, n, size=n)
#生成0-1000的随机整数，作为数据序列位置号。
#此函数精髓就是利用数组位置号来再抽样
data_sample = array[index_arr]
#根据生成的随机数据序列号，来再次抽样赋值新数组。
sample_result = func(data_sample)
sample_result_arr.append(sample_result)
#append函数是将重采样的sample_result数据添加到sample_result_arr数组后面
a = 1 - c
#如果是95%置信度，那就是 c = 0.95 ,a = 1- c= 0.05
k1 = int(B * a / 2)
#如果B是1000，则k1代表从小到大排列，第2.5%个分位处的序列位置号。
k2 = int(B * (1 - a / 2))
#如果B是1000，则k1代表从小到大排列，第97.5%个分位处的序列位置号。
auc_sample_arr_sorted = sorted(sample_result_arr)
#将1000次重新抽样的数据从大到小排列。
lower = auc_sample_arr_sorted[k1]
#取2.5%分位处的值
higher = auc_sample_arr_sorted[k2]
#取2.5%分位处的值
return lower, higher
#返回上下置信度。
#-----------自定义函数结束
if __name__ == '__main__':
result = bootstrap(np.random.randint(0, 50, 50), 1000, 0.95, average)
#构建0-50之间的50个随机数，1000次重采样，95%置信度，平均值输出
print(result)
#输出结果；平均值信度区间为(CI:20.48, 28.32)
```**

**

NCL代码实现：

**

begin
data   = random_uniform(0, 50, 50)  ;创建0-50的50个随机数据
;======================================================================
;计算置信度区间
;======================================================================
z     = data   ;原始数据
nBoot = 1000   ;重复提取次数
c     = 0.95   ;置信度
nDim  = 0      ;维度
opt   = False
start = bootstrap_optarg(z, nBoot, nDim, opt)   ;提取所需数据
dimz    = start[0]
rankz   = start[1]
N       = start[2]
NN      = start[3]      ; NN=N (most commonly)
zBoot   = start[4]
delete(start)
do nb=0,nBoot-1
iw = bootstrap_indices(N,NN,opt)                ;重新生成随机序列号
newdata = z(iw)                                 ;重建数组
zBoot(nb) = dim_avg_n_Wrap(newdata, 0)   ;利用重建数据求平均，一共求1000次
end do
c1 = 1 - c
;#如果是95%置信度，那就是 c = 0.95 ,a = 1- c= 0.05
k1 = floattoint(nBoot * c1 / 2)
;#如果nBoot是1000，则k1代表从小到大排列，第2.5%个分位处的序列位置号。
k2 = floattoint(nBoot * (1 - c1 / 2))
;#如果nBoot是1000，则k1代表从小到大排列，第97.5%个分位处的序列位置号
qsort(zBoot)
;#将1000次重新抽样的数据从大到小排列。
lower = zBoot(k1)                                 ;#取2.5%分位处的值
higher = zBoot(k2)                                ;#取97.5%分位处的值
print(lower)                                      ;输出到屏幕
print(higher)
;===========================================================================================
end
;=========================================================================

;---------------------------------

2. GEV拟合PDF的bootstrap

NCL代码实现：

begin
f1       = "G:/attribution/cmip5/EAWM.txt"  ;读取数据
data     = asciiread(f1,(/ 525/),"float")
qsort(data)  ;从小到大排序
x2018    = -43.16885          ;2018年当年的值， 从图片读出来的概率为0.57%
;---GEV parameter estimation
;****************************************** 计算原始数据中指定阈值的GEV-pdf值
N        =  100
xMin     = min(data )
xMax     = max(data )
x1       = fspan(xMin,xMax,N)
vals     = extval_mlegev (data  , 0, False)    ; dims=0
center   = vals(0)                            ; extract for clarity
scale    = vals(1)
shape    = vals(2)
pdfcdf   = extval_gev(x1 , shape, scale, center,0)
pdf      = pdfcdf[0]                          ; extract from list for convenience
delete(pdfcdf)
if (scale.le.0)  then
print("extval_gev: unexpected value(s): scale="+scale+" shape="+shape)
continue
end if
y   = (x2018-center)/scale
a   = 1+shape*y
q   = 1.0/scale
if (shape.ne.0) then
p1      = -1-(1.0/shape)
p2      = -1.0/shape
pdf2018 = q*exp(-(a^p2))*(a^p1)
else
pdf2018 = q*exp(-y-exp(-y))
pdf2018 = where(pdf2018.lt.0, pdf2018@_FillValue, pdf2018)
end if
print(pdf2018)                       ;========原始数据中指定阈值的GEV-pdf值
;***************************************************************

;======================================================================
;以下是指定阈值的GEV-pdf的bootstrap,计算置信度区间 0.57% (CI:0.47%-0.63%)
;======================================================================
z     = data   ;原始数据
z2018 = x2018  ;当年阈值
nBoot = 1000   ;重复提取次数
c     = 0.95   ;置信度
nDim  = 0      ;维度
opt   = False

start = bootstrap_optarg(z, nBoot, nDim, opt)   ;提取所需数据
dimz    = start[0]
rankz   = start[1]
N       = start[2]
NN      = start[3]      ; NN=N (most commonly)
zBoot   = start[4]
delete(start)

do nb=0,nBoot-1
iw = bootstrap_indices(N,NN,opt)                ;重新生成随机序列号
newdata = z(iw)                                 ;重建数组
;--------------------------------------------------------------------------计算指定阈值的GEV值
vals     = extval_mlegev (newdata  , 0, False)
center   = vals(0)
scale    = vals(1)
shape    = vals(2)
delete(vals)
delete(newdata)
if (scale.le.0)  then
print("extval_gev: unexpected value(s): scale="+scale+" shape="+shape)
continue
end if
y   = (z2018-center)/scale
a   = 1+shape*y
q   = 1.0/scale
if (shape.ne.0) then
p1       = -1-(1.0/shape)
p2       = -1.0/shape
pdf2018  = q*exp(-(a^p2))*(a^p1)
else
pdf2018  = q*exp(-y-exp(-y))
pdf2018  = where(pdf2018.lt.0, pdf2018@_FillValue, pdf2018)
end if
delete(center)
delete(scale)
delete(shape)
delete(y)
delete(a)
delete(q)
delete(p1)
delete(p2)
zBoot(nb) = pdf2018                            ;指定阈值对应的PDF值，将重复计算的对应阈值的PDF存储到数组zBoot
delete(pdf2018)
end do

c1 = 1 - c
;#如果是95%置信度，那就是 c = 0.95 ,a = 1- c= 0.05
k1 = floattoint(nBoot * c1 / 2)
;#如果nBoot是1000，则k1代表从小到大排列，第2.5%个分位处的序列位置号。
k2 = floattoint(nBoot * (1 - c1 / 2))
;#如果nBoot是1000，则k1代表从小到大排列，第97.5%个分位处的序列位置号。
qsort(zBoot)
;#将1000次重新抽样的数据从大到小排列。
lower = zBoot(k1)
;#取2.5%分位处的值
higher = zBoot(k2)
;#取97.5%分位处的值
print(lower)
;输出到屏幕
print(higher)
;===========================================================================================
;===========================================================================================
end

GEV拟合pdf的Python代码实现：

正在加载中。。。。。。。。。。。

References:

1. 统计中置信区间介绍
2. Bootstrap方法介绍
3. NCL官网：GEV函数
4. NCL官网：Bootstrap相关函数
5. 极大似然估计原理1（即GEV拟合原理）
6. 极大似然估计原理2
7. 极大似然估计python实现
8. GEV拟合python代码实现