如何在一个图中绘制拟合图和伽马分布的实际图？

Posted 2023-03-12

【中文标题】如何在一个图中绘制拟合图和伽马分布的实际图？

【英文标题】：How to draw fitted graph and actual graph of gamma distribution in one plot?

【发布时间】：2016-04-16 17:29:04

【问题描述】：

步骤 1。加载所需的包。

library(ggplot2)
library(MASS)

步骤 2。生成 10,000 个适合 gamma 分布的数字。

x <- round(rgamma(100000,shape = 2,rate = 0.2),1)
x <- x[which(x>0)]

步骤 3。画出pdf（概率密度函数），假设我们不知道x适合哪个分布。

t1 <- as.data.frame(table(x))
names(t1) <- c("x","y")
t1 <- transform(t1,x=as.numeric(as.character(x)))
t1$y <- t1$y/sum(t1[,2])
ggplot() + geom_point(data = t1,aes(x = x,y = y)) + theme_classic()

步骤 4。从图中我们可以看出x的分布很像gamma分布，所以我们使用MASS包中的fitdistr()来获取gamma分布的shape和rate的参数。

fitdistr(x,"gamma") 
##       output 
##       shape           rate    
##   2.0108224880   0.2011198260 
##  (0.0083543575) (0.0009483429)

第 5 步。在同一个图中画出实际点（黑点）和拟合图（红线），问题来了，请先看图。

ggplot() + geom_point(data = t1,aes(x = x,y = y)) +     
geom_line(aes(x=t1[,1],y=dgamma(t1[,1],2,0.2)),color="red") + 
theme_classic()

问题一：真正的参数是shape=2，rate=0.2，我用函数fitdistr()得到的参数是shape=2.01，rate=0.20。这两个几乎一样，但是为什么拟合图不能很好地拟合实际点，肯定是拟合图有问题，或者我绘制拟合图和实际点的方式完全错误，我该怎么办?

问题2：得到我建立的模型的参数后，我以哪种方式评估模型，例如RSS(residual square sum) 用于线性模型，或者p-value 的shapiro.test() ，ks.test() 和其他测试?我统计知识贫乏，请大家帮我解答一下，谢谢！（ps：我在google和***上搜索过很多次，但都没有用，所以不要投这个问题没用，谢谢！ )

【问题讨论】：

我没有太多的声誉可以直接显示图片，有人可以帮我一个忙吗？感谢您的帮助！

请显示您在第 5 步中使用的代码（即绘制数据和曲线）。

对不起，我忘记了，请再次查看描述@jbaums

你提出的问题实际上不在Stack Overflow的范围内，真正属于Cross Validated。对于您的 Q1，我认为您将离散箱的相对频率（概率质量函数）与连续概率密度混淆：有关一些讨论，请参阅 here、here 和 here。关于 Q2，这不是统计问题的地方（如果您有关于实施特定程序的编程问题，可以）。试试Cross Validated的人。

@jbaums 我已经解决了我关于简历的问题，谢谢你的帮助

【参考方案1】：

由于某种原因，您的拟合线似乎恰好高了 10 倍：

ggplot() + geom_point(data = t1,aes(x = x,y = y)) +     
  geom_line(aes(x=t1[,1],y=(dgamma(t1[,1],2,0.2))/10),color="red") + 
  theme_classic()

完美契合：

正如 jbaums 所说，这是由每个 dx 的密度在这种情况下为 0.1 引起的。

【讨论】：

密度以每 dx 的密度给出，在这种情况下 dx 为 0.1，因为 OP 将数据四舍五入到小数点后 1 位。一般情况下，density = count/(sum(count)*binwidth).

@RHA 是啊，这个方法我已经试过了，但是就像你说的，为什么？

@jbaums 能不能重写代码，你说的我完全听不懂，谢谢

@jbaums 数据主要用于学术研究，所以我想知道内在原因，不只是为了回答问题，谢谢

@LingZhang - 查看我对您帖子的评论。