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.05k1 = 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-1iw = 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) thenprint("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) thenp1 = -1-(1.0/shape) p2 = -1.0/shapepdf2018 = q*exp(-(a^p2))*(a^p1) elsepdf2018 = 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 = Falsestart = 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-1iw = 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) thenprint("extval_gev: unexpected value(s): scale="+scale+" shape="+shape)continueend ify = (z2018-center)/scale a = 1+shape*yq = 1.0/scaleif (shape.ne.0) thenp1 = -1-(1.0/shape) p2 = -1.0/shapepdf2018 = q*exp(-(a^p2))*(a^p1) elsepdf2018 = q*exp(-y-exp(-y))pdf2018 = where(pdf2018.lt.0, pdf2018@_FillValue, pdf2018)end ifdelete(center)delete(scale)delete(shape)delete(y)delete(a)delete(q)delete(p1)delete(p2)zBoot(nb) = pdf2018 ;指定阈值对应的PDF值,将重复计算的对应阈值的PDF存储到数组zBootdelete(pdf2018)
end doc1 = 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:
- 统计中置信区间介绍
- Bootstrap方法介绍
- NCL官网:GEV函数
- NCL官网:Bootstrap相关函数
- 极大似然估计原理1(即GEV拟合原理)
- 极大似然估计原理2
- 极大似然估计python实现
- GEV拟合python代码实现
Bootstrap置信区间和GEV拟合pdf相关推荐
- bootstrap table export插件导出pdf格式文件中文乱码问题解决办法
bootstrap table export插件导出pdf格式文件中文乱码的问题折腾了我整整两天,网上到处都是改源码,自己设置字体的方案,我都没搞定.结果今天看到官方文档(地址:GitHub - hh ...
- php 置信区间 计算,bootstrap置信区间如何求
bootstrap置信区间: 假设总体的分布F未知,但有一个容量为n的来自分布F的数据样本,自这一样本按有放回抽样的方法抽取一个容量为n的样本,这种样本称为bootstrap样本.相继地.独立地自原始 ...
- matlab 提取极值,利用matlab 进行极值统计的一个例子——gev 方法.pdf
利用matlab 进行极值统计的一个例子--gev 方法 利用 Matlab 进行极值统计的一个例子--GEV 方法 科研菜鸟 /u/sanshiphy 数据和例子均来自于 S. Coles, An ...
- matlab对闭合轮廓进行多边形逼近,物体轮廓线的多边形拟合.PDF
物体轮廓线的多边形拟合 物体轮廓线的多边形拟合 1 2 魏炎初 张建军 (清华大学信息网络工程研究中心) 摘要:本文讨论了在给定数字轮廓线顶点数目为 N 的情况下,从其上选择 k 个点来构造拟和 多边 ...
- python广义极值_广义极值(GEV)极大似然拟合的奇异pdf
首先,我认为您可能希望将位置参数固定在0.在 其次,你的数据中有0,结果拟合可能在x=0处有{}pdf,例如对于GEV拟合或Weibull拟合. 因此,拟合实际上是正确的,但是当绘制pdf(包括x=0 ...
- bootstrap使用tableExport导出pdf时中文乱码问题
前言 最近拿到了一个任务,让处理一下公司的系统平台问题.问题就是页面导出PDF文件,有中文的话显示的都是乱码.因为公司的项目都是给国外客户使用的,所以我估计从设计到测试都没有考虑中文的问题.但是为啥现 ...
- Bootstrap重采样进行参数估计 - 置信区间
Bootstrap重采样进行参数估计 - 置信区间 文章目录 Bootstrap重采样进行参数估计 - 置信区间 一.Bootstrap简介 二.为什么要使用Bootstrap 三.经验Bootstr ...
- Cox 比例风险模型中HR和置信区间
Cox 回归是一种用于生存分析的统计模型,它可以用来估计某个因素对生存时间的影响.Cox 回归基于 Cox 比例风险模型,该模型假设风险比率是常数,即不随时间变化.在 Cox 回归中,我们使用最大似然 ...
- 深度解析机器学习中的置信区间(附代码)
作者:Jason Brownlee 翻译:和中华 校对:丁楠雅 本文约4000字,建议阅读15分钟. 本文介绍了置信区间的概念以及如何计算置信区间和bootstrap置信区间. 机器学习很多时候需要估 ...
最新文章
- 表格对决CSS--一场生死之战
- 转:SLAM算法解析:抓住视觉SLAM难点,了解技术发展大趋势
- qt关于添加模块的说明
- JZOJ 2413. 【NOI2005】维护数列
- 《SAS编程与数据挖掘商业案例》学习笔记之四
- [js] 使用delete删除数组,其长度会改变吗
- gitlab ci 配置 java_GitLab CI/CD 配置
- python c语言接口_C/C++ 提供 Python 接口
- OpenStack源码系列---nova-conductor
- [转]PKM2:优秀的个人知识管理工具
- 查看与设置华为路由器的版本信息,名称、时钟、登录提示
- C++STL算法equal(15)
- Makefile wildcard函数说明
- Unity3D 参数曲线 实现曲线上的匀速运动
- win7、win10系统双屏显示任务栏
- java图书管理系统这个怎么改呢
- CSS3小可爱亲吻表白特效,给你的五一假期增添点小乐趣
- Rosalind Java|Finding a Spliced Motif
- HTTP协议与Get和Post的区别
- 项目管理的10大知识领域之范围管理