如何使用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)

我不确定我是否以正确的方式进行操作。任何指导都将不胜感激。

- Chris

---

你从哪里获取了 dgumbel？它不是 r-base 的分发版。即使已经加载了 VGAM，gumbel() 调用仍然会抛出一个错误。 - IRTFM

@DWin 从包evd中获取dgumbel（参见：http://tiny.cc/8izrk）。是否有基本函数？ - Chris

在evd包的内容中没有“gumbel”函数。 - IRTFM

@DWin gumbel是evir包的一个函数。但问题实际上并不是关于这个函数，因为我正在尝试绘制属于evd包的函数dgumbel，对吧？ - Chris

如果问题与evir函数覆盖evd函数无关，那么这可能是正确的。尝试进行干净的会话。不要加载evir（它没有命名空间，因此您无法使用:::运算符），然后运行您的绘图。另外，您应该编辑您的问题以指示需要哪些库或require调用。 - IRTFM

应该说是evd没有命名空间。 - IRTFM

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'))

- Ramnath

---

+1 很棒的答案。自ggplot2版本0.9.1起，“ops”已被弃用，请改用“theme”：theme（legend.position ='top'，legend.direction ='horizontal'） - Eduardo

Answer 2

之前的会话显示，来自gumbel call的参数估计值接近于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可以避免与evir中的dgumbel相关函数发生冲突。当您第二次加载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))

我认为添加经验“密度”的建议是不错的，但我更喜欢使用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))

- IRTFM

---

好的，你说的包冲突是正确的。我已经编辑了问题以反映所需的包。我希望将我的数据拟合到gumbel分布，因此尝试使用evir中的函数。这有意义吗？有没有其他选择？ - Chris

我理解。事实上，我利用之前会话中使用“gumbel”的知识，替换了更有意义的值来调用“dgumbel”。你可能需要运行“gumbel”调用，注意参数估计值，然后分离evir包，加载evd包并进行绘图。或者您可以使用具有NAMESPACE的dgumbel软件包。 - IRTFM

我尝试了ismeve包中的gum.fit，它似乎也能胜任这项工作。感谢您解决了这个问题！ - Chris

@DWin。你的第一个解决方案没有绘制经验密度图。你可能需要这样调用它：ggplot(data=speeds2, aes(x=speed)) + geom_density() + stat_function(fun=dgumbel, args=list(loc=24, scale=11))。 - Ramnath

@Ramnath。同意。添加了代码来做类似的事情，但它的平滑程度不如geom_density的默认设置。 - IRTFM

Answer 3

通过对代码进行小的修改（添加一个几何对象），我成功地解决了问题。

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

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

- Richie Cotton