0

我想在观察数据的直方图上绘制从一组观察结果中得出的伽马密度函数。我能够为伽马拟合生成直方图和参数估计。这是针对来自主数据集的多个数据子集完成的。如何在此循环中创建的每个直方图上绘制伽马密度函数？

我目前有：

library(MASS)

species <- c('acesac', 'acesac', 'acesac', 'acesac', 'acesac', 'acesac',
 'polbif', 'polbif', 'polbif', 'polbif', 'polbif', 'polbif')
tmean <- c(2,3,5,6,6,7,5,6,6,6,8,9) 
Data <- data.frame(species, tmean) 

for (i in unique(Data$species)){
  subdata <- subset(Data, species ==i)
  hist(subdata$tmean, main = i)
  dist <- fitdistr(subdata$tmean, "gamma")
}

我在想我应该使用lines()，但是，不确定如何指定这个？

4

1 回答 1

0

我会添加library(MASS)到您的示例中。您可能想尝试使用curve和使用add = TRUE. 另一种选择是使用它，library(fitdistrplus)因为它可以直接绘制输出dist；但是，我找不到（快速）更改情节标题的方法。

library(MASS)
species <- c('acesac', 'acesac', 'acesac', 'acesac', 'acesac', 'acesac',
             'polbif', 'polbif', 'polbif', 'polbif', 'polbif', 'polbif')
tmean <- c(2,3,5,6,6,7,5,6,6,6,8,9) 
Data <- data.frame(species, tmean) 

for (i in unique(Data$species)) {
  subdata <- subset(Data, species ==i)
  dist <- fitdistr(subdata$tmean, "gamma")
  hist(subdata$tmean, main = i)
  curve(dgamma(x, shape = dist$estimate[1], rate = dist$estimate[2]), 
        add = TRUE,
        col = "red")
}

根据我关于library(fitdistrplus)查看以下输出的评论：

library(fitdistrplus)

for (i in unique(Data$species)) {
  subdata <- subset(Data, species ==i)
  dist <- fitdistrplus::fitdist(subdata$tmean, "gamma")
  plot(dist)
}

请注意，您会获得其他图表（QQ 图、经验和理论 CFS 以及 PP 图），并且“免费”绘制了密度线。但是，您将失去添加main = i. 我敢肯定，比我聪明的人可以想出一种快速添加标题或修改方法的fitdistrplus::plot.fitdist方法-可能值得提出一个单独的问题。

于 2015-10-09T19:21:44.423 回答