将分布拟合到经验数据

Question

我正在尝试将 beta 分布拟合到根据经验数据创建的直方图。

我遇到的问题是拟合分布远高于原始直方图中的条形。

原始数据超出 [0,1] 的范围，这是可以评估 beta 分布的范围，因此我重新调整原始数据，使其位于区间 [0,1] 内。

这是我的代码：

 load("https://www.dropbox.com/s/c3psxx8jjbc20mo/data.Rdata?dl=0")

  #create histogram with values normalized between 0 and 1
  h <- hist((data-min(data)) / (max(data)-min(data)),lty="blank",col="grey")
  #normalize the density so the y-axis goes from 0 to 1
  h$density <- h$counts/max(h$counts)
  #plot the results
  plot(h,freq=FALSE,cex.main=1,cex.axis=1,yaxt='n',ylim=c(0,1.5),col='grey',lty='blank',xaxt='n')
  axis(2,at=seq(0,1,0.5),labels=seq(0,1,0.5))
  axis(1,at=seq(0,1,0.5),labels=seq(0,1,0.5))

  #fit beta distribution
  a <- (data-min(data)) / (max(data)-min(data))
  a[a==1] <- 0.9999
  a[a==0] <- 0.0001
  fit.beta <- suppressWarnings(fitdistr(a, "beta", start = list( shape1=0.1, shape2=0.1 ) ))

  #overlay curve from beta distribution
  alpha <- fit.beta$estimate[1]
  beta <- fit.beta$estimate[2]
  b <- rbeta(length(data),alpha,beta)
  lines(density(b))

我错过了什么？

Answer 1

首先，您需要使用hist(..., freq=TRUE) 作为直方图。然后，要正确设置 y 轴范围，您可以计算 beta 分布的最大值 (see e.g. here)。最后，使用dbeta 比生成随机样本然后估计密度要好得多：

maxibeta <- dbeta((alpha-1)/(alpha+beta-2), alpha, beta)
hist( (data-min(data)) / (max(data)-min(data)), 
      prob=TRUE, col="grey", border="white", ylim=c(0, maxibeta), 
      main="Histogram + fitted distribution")
plot(function(x) dbeta(x,alpha,beta), add=TRUE, col=2, lwd=2)

---

编辑：一个更通用的解决方案，但这让我有点难过，因为它没有使用 beta 发行版的好特性：

fbeta <- function(x)  dbeta(x,alpha,beta)
maxibeta <- optimize(fbeta, interval = c(0,1), maximum = TRUE)$objective

histo <- hist((data-min(data)) / (max(data)-min(data)), plot = FALSE)

plot(histo, freq=FALSE, col="grey", border="white", 
     ylim=c(0, max(maxibeta, max(histo$density))), 
     main="Histogram + fitted distribution")
plot(fbeta, add=TRUE, col=2, lwd=2)