2013-03-26 111 views
2

我有風數據用於執行極值分析(計算回報水平)。我正在使用R,套件包括'evd','extRemes'和'ismev'。威布爾分佈的置信區間

我裝修GEV,岡貝爾和Weibull分佈,以便估計回報水平(RL)的一些週期T
對於極值分佈和岡貝爾情況下,我可以使用極端得到RL的和置信區間:: return.level()函數。

一些代碼:

require(ismev) 
require(MASS) 

data(wind) 
x = wind[, 2] 
rperiod = 10 

fit <- fitdistr(x, 'weibull') 
s <- fit$estimate['shape'] 
b <- fit$estimate['scale'] 

rlevel <- qweibull(1 - 1/rperiod, shape = s, scale = b) 

## CI around rlevel 
## ci.rlevel = ?? 

不過對於韋伯的情況下,我需要一些幫助,從而生成CI的。

+1

你已經提供了一個粗略的描述你做了什麼,但這是一個編碼論壇,你需要發佈數據和代碼。 – 2013-03-26 15:46:23

+0

添加代碼,謝謝。 – Fernando 2013-03-26 15:53:52

回答

2

我懷疑令人難以忍受的正確答案是聯合置信區域是一個橢圓或一些彎曲的香腸形狀,但是您可以使用vcov函數從擬合對象中提取參數的方差估計,然後構建標準誤+/- 1.96 SE的應翔實:

> sqrt(vcov(fit)["shape", "shape"]) 
[1] 0.691422 
> sqrt(vcov(fit)["scale", "scale"]) 
[1] 1.371256 

> s +c(-1,1)*sqrt(vcov(fit)["shape", "shape"]) 
[1] 6.162104 7.544948 
> b +c(-1,1)*sqrt(vcov(fit)["scale", "scale"]) 
[1] 54.46597 57.20848 

的常用方法來計算一個參數一個CI是假設正態分佈和使用THETA +/- 1.96 * SE(THETA)。在這種情況下,你有兩個參數,兩個參數都會給你一個「盒子」,一個間隔的2D模擬。真正正確的答案在'scale'by-'hape'參數空間中會更復雜一些,並且可能最容易用模擬方法實現,除非您比我更好地理解理論。

+0

所以現在我可以使用這些值來計算qweibull_inf和qweibull_hi? – Fernando 2013-03-28 16:56:56

+0

我不完全確定我在這裏理解你的術語。 – 2013-03-28 17:04:13

+0

我有形狀和規模的CI。我可以用qweibull(p,shape,scale)來計算weibull分位數嗎? – Fernando 2013-03-28 17:24:59