1

我必须将伽玛分布曲线叠加到其他类似幂律曲线的图中。我首先以对数刻度绘制直方图的点

  plot(log(pp$mids),log(pp$density))

然后我想叠加我的伽玛分布曲线调用一个外部函数 gamma()

  gamma <- function(X)
  {
  n <- length(X)
  theta<-var(hh2$V1)/mean(hh2$V1)
  kappa<-mean(hh2$V1)/theta
  y<-rgamma(n,kappa,theta)
  xx<-hist(y,plot=F)
  curve(log(xx$density),add=T,col='violet',type='l')
  return( c(kappa) ) 
  }

但这会给我一个错误，因为 curve() 需要一条真正的曲线来绘制。我怎样才能做到这一点？

4

1 回答 1

1

这是您的代码的一个有点工作的变体：

生成一个结构为（我猜）您的数据的示例：

library(rmutil)  ## for rpareto
set.seed(101)
hh2 <- data.frame(V1=rpareto(1000, m=1, s=1.5))

初始直方图计算：

pp <- hist(hh2$V1,plot=FALSE)

函数（最好不要调用它gamma，因为它掩盖了内置函数）：

ghistfun <- function(x) {
    n <- length(x)
    scalepar <- var(x)/mean(x)
    shapepar <- mean(x)^2/var(x)
    y <- rgamma(n,shape=shapepar,scale=scalepar)
    xx <- hist(y,plot=FALSE)
    lines(log(xx$mids),log(xx$density),col="red")
    curve(dgamma(exp(x),shape=shapepar,scale=scalepar,log=TRUE),
        add=TRUE,col="blue")

    shapepar
}

使用一个非常大的数字而不是仅仅使用数据的长度可能会更好n，除非您对查看完全相同大小的数据集中的随机波动特别感兴趣。或者，您可以只使用curve(dgamma(x,...))，如图所示（我最初认为您必须允许将密度从 缩放x到log(x)，但是由于您在原始（未记录）比例上计算直方图然后转换的方式bin 中点，您不必...）

plot(log(pp$mids),log(pp$density))
ghistfun(hh2$V1)

在此处输入图像描述

于 2012-07-27T15:06:11.087 回答