R:如何使用ggplot2的stat_function绘制gumbel分布

Question

如果这是相当脆弱的话,请耐心等待,如果我遗漏了任何东西,请随时提问...

我试图根据以下链接进行50年的极端风计算

http://www.wasp.dk/Products/weng/ExtremeWinds.htm

他们似乎使用gumbel分布,所以我在包"evir"中使用了函数gumbel以适应数据的分布,并在包"evd"中使用dgumbel作为绘图函数.

package("evd")
package("evir")

speeds2 <- data.frame(speed=sample(10:50,1000,rep=TRUE))
gumbel(speeds2$speed)

然后我尝试使用ggplot2的stat_function来绘制它,就像这样(除了现在我已经为loc和scale添加了虚拟值).

library(ggplot2)
ggplot(data=speeds2, aes(x=speed)) + 
  stat_function(fun=dgumbel, args=list(loc=1, scale=0.5))

我收到以下错误:

Error in dgev(x, loc = loc, scale = scale, shape = 0, log = log) : 
  unused argument(s) (loc = loc, scale = scale, shape = 0, log = log)

我不确定我是否以正确的方式这样做.任何指针都将非常感激.

Answer 1

这是我编写的一个通用函数,用于简化具有拟合和经验密度的数据绘图.

# FUNCTION TO DRAW HISTOGRAM OF DATA WITH EMPIRICAL AND FITTED DENSITITES
# data  = values to be fitted
# func  = name of function to fit (e.g., 'norm', 'gumbel' etc.)
# start = named list of parameters to pass to fitting function 
hist_with_density = function(data, func, start = NULL){
    # load libraries
    library(VGAM); library(fitdistrplus); library(ggplot2)

    # fit density to data
    fit   = fitdist(data, func, start = start)
    args  = as.list(fit$estimate)
    dfunc = match.fun(paste('d', func, sep = ''))

    # plot histogram, empirical and fitted densities
    p0 = qplot(data, geom = 'blank') +
       geom_line(aes(y = ..density..,colour = 'Empirical'),stat = 'density') +
       stat_function(fun = dfunc, args = args, aes(colour = func))  +
       geom_histogram(aes(y = ..density..), alpha = 0.4) +
       scale_colour_manual(name = '', values = c('red', 'blue')) + 
       opts(legend.position = 'top', legend.direction = 'horizontal')
    return(p0)  
}

以下是如何使用它的两个示例示例1:适合Gumbel

data1 = sample(10:50,1000,rep=TRUE)
(hist_with_density(data1, 'gumbel', start = list(location = 0, scale = 1)))

示例2:拟合正态分布

data2 = rnorm(1000, 2, 1)
(hist_with_density(data2, 'norm'))

Answer 2

早些时候的会议显示，gumbel 调用的参数估计值接近 24 和 11。

library(evd)
library(ggplot2)
 speeds2 <- data.frame(speed=sample(10:50,1000,rep=TRUE))
 ggplot(data=speeds2, aes(x=speed), geom="density") + 
   stat_function(fun=dgumbel, args=list(loc=24, scale=11))

如果仅使用 1 和 0.5 的参数，您会得到一条平直的线。加载仅evd防止与 中 dgumbel 相关函数发生冲突evir。当您加载evir第二个时，您会得到：

> speeds2 <- data.frame(speed=sample(10:50,1000,rep=TRUE))
> ggplot(data=speeds2, aes(x=speed), geom="density") + 
+   stat_function(fun=dgumbel, args=list(loc=24, scale=11))
Error in dgev(x, loc = loc, scale = scale, shape = 0, log = log) : 
  unused argument(s) (loc = loc, scale = scale, shape = 0, log = log)

演示如何调用dgumbel特定（表现更好）包中的函数：

library(VGAM)
ggplot(data = speeds2, aes(x = speed)) + 
   stat_function(fun = VGAM::dgumbel, args = list(location = 24, scale = 11))

我认为 Ramnath 添加经验“密度”的建议很好，但我更喜欢使用 geom_histogram：

ggplot(data=speeds2, aes(x=speed)) + geom_histogram(aes(y = ..density..) , binwidth=5 ) + 
                            stat_function(fun=dgumbel, args=list(loc=24, scale=11))