R：如何使用 ggplot2 的 stat_function 绘制 gumbel 分布答案

【问题标题】：R: How to plot gumbel distribution using ggplot2's stat_functionR：如何使用 ggplot2 的 stat_function 绘制 gumbel 分布
【发布时间】：2011-10-14 10:00:03
【问题描述】：

如果这很脆弱，请多多包涵，如果我遗漏了什么，请随时提出问题...

我正在尝试根据以下链接进行一些 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)

我不确定我这样做是否正确。任何指针将不胜感激。

【问题讨论】：

你从哪里得到dgumbel？它不是基于 r 的发行版。即使加载了 VGAM，gumbel() 调用也会引发错误。
@DWin 从包evd 中获取dgumbel（参见：tiny.cc/8izrk）。有基础函数吗？
在 package:evd 的内容中没有 gumbel 函数。
@DWin gumbel 是包 evir 的一个功能。问题不在于这个函数，因为我实际上是在尝试绘制函数dgumbel，它确实属于包evd，不是吗？
如果问题与覆盖evd 函数的evir 函数无关，则可能是这样。尝试干净的会话。不要加载 evir （它没有 NAMESPACE 所以你不能使用 ':::' 操作符。然后运行你的情节。你也应该编辑你的问题以表明需要什么库或需要调用。

标签： r ggplot2

【解决方案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'))

【讨论】：

+1 很好的答案。 “ops”自 ggplot2 0.9.1 版起已弃用，请改用“主题”：theme(legend.position = 'top', legend.direction = 'horizontal')

【解决方案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 可防止与 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))

我认为 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))

【讨论】：

好的，关于包冲突你是对的。我已经编辑了问题以反映所需的包。我希望将 gumbel 分布拟合到我的数据中，因此尝试使用 evir 中的函数。那有意义吗？有其他选择吗？
对我来说很有意义。事实上，我使用从早期会话使用gumbel 中获得的知识来替换dgumbel 调用的更有意义的值。您可能需要运行“gumbel”调用，注意参数估计，然后分离包 evir，然后加载包 evd 并进行绘图。或者您可以使用具有 NAMESPACE 的 dgumbel 包。
我在 ismev 包中尝试了 gum.fit，它似乎也可以完成这项工作。感谢您解决这个问题！
@DWin。您的第一个解决方案没有绘制经验密度。您可能需要将其称为ggplot(data=speeds2, aes(x=speed)) + geom_density() + stat_function(fun=dgumbel, args=list(loc=24, scale=11))。
@Ramnath。同意。添加了代码来做类似的事情，除了它不像 geom_density 的默认设置那样平滑。

【解决方案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()

【讨论】：