Comment dessiner un graphique ajusté et un graphique réel de la distribution gamma dans un seul graphique?

Chargez le package nécessaire.

library(ggplot2)
library(MASS)

Générez 10 000 nombres adaptés à la distribution gamma.

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

Dessinez la fonction de densité de probabilité, supposant que nous ne savons pas à quelle distribution x correspond.

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()

À partir du graphique, nous pouvons apprendre que la distribution de x est assez similaire à la distribution gamma, nous utilisons donc fitdistr()dans le package MASSpour obtenir les paramètres de forme et de taux de distribution gamma.

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

Dessinez le point réel (point noir) et le graphique ajusté (ligne rouge) dans le même tracé, et voici la question, veuillez d'abord regarder le tracé.

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()

J'ai deux questions:

1. Les paramètres réels sont shape=2, rate=0.2et les paramètres que j'utilise la fonction fitdistr()pour obtenir sont shape=2.01, rate=0.20. Ces deux sont presque les mêmes, mais pourquoi le graphique ajusté ne correspond pas bien au point réel, il doit y avoir quelque chose de mal dans le graphique ajusté, ou la façon dont je dessine le graphique ajusté et les points réels est totalement erronée, que dois-je faire ?

2. Après avoir obtenu le paramètre du modèle que j'établis, de quelle manière j'évalue le modèle, quelque chose comme RSS (somme carrée résiduelle) pour le modèle linéaire, ou la valeur de p de shapiro.test(), ks.test()et un autre test?

Je suis pauvre en connaissances statistiques, pourriez-vous bien vouloir m'aider?

ps: J'ai fait des recherches dans Google, stackoverflow et CV plusieurs fois, mais je n'ai rien trouvé en rapport avec ce problème

question 1

La façon dont vous calculez la densité à la main semble incorrecte. Il n'est pas nécessaire d'arrondir les nombres aléatoires de la distribution gamma. Comme l'a noté @Pascal, vous pouvez utiliser un histogramme pour tracer la densité des points. Dans l'exemple ci-dessous, j'utilise la fonction densitypour estimer la densité et la représenter sous forme de points. Je présente l'ajustement à la fois avec les points et avec l'histogramme:

library(ggplot2)
library(MASS)

# Generate gamma rvs

x <- rgamma(100000, shape = 2, rate = 0.2)

den <- density(x)

dat <- data.frame(x = den$x, y = den$y)

# Plot density as points

ggplot(data = dat, aes(x = x, y = y)) + 
  geom_point(size = 3) +
  theme_classic()

# Fit parameters (to avoid errors, set lower bounds to zero)

fit.params <- fitdistr(x, "gamma", lower = c(0, 0))

# Plot using density points

ggplot(data = dat, aes(x = x,y = y)) + 
  geom_point(size = 3) +     
  geom_line(aes(x=dat$x, y=dgamma(dat$x,fit.params$estimate["shape"], fit.params$estimate["rate"])), color="red", size = 1) + 
  theme_classic()

# Plot using histograms

ggplot(data = dat) +
  geom_histogram(data = as.data.frame(x), aes(x=x, y=..density..)) +
  geom_line(aes(x=dat$x, y=dgamma(dat$x,fit.params$estimate["shape"], fit.params$estimate["rate"])), color="red", size = 1) + 
  theme_classic()

Voici la solution fournie par @Pascal:

h <- hist(x, 1000, plot = FALSE)
t1 <- data.frame(x = h$mids, y = h$density)

ggplot(data = t1, aes(x = x, y = y)) + 
  geom_point(size = 3) +     
  geom_line(aes(x=t1$x, y=dgamma(t1$x,fit.params$estimate["shape"], fit.params$estimate["rate"])), color="red", size = 1) + 
  theme_classic()

question 2

Pour évaluer la qualité de l'ajustement, je recommande le package fitdistrplus. Voici comment il peut être utilisé pour ajuster deux distributions et comparer leurs ajustements graphiquement et numériquement. La commande gofstataffiche plusieurs mesures, telles que l'AIC, le BIC et certaines statistiques de gof telles que le test KS, etc. Elles sont principalement utilisées pour comparer les ajustements de différentes distributions (dans ce cas, gamma par rapport à Weibull). Plus d'informations peuvent être trouvées dans ma réponse ici :

library(fitdistrplus)

x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00,
       38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40,
       42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40,
       49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60,
       45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30,
       36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00,
       38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34)

fit.weibull <- fitdist(x, "weibull")
fit.gamma <- fitdist(x, "gamma", lower = c(0, 0))

# Compare fits 

graphically

par(mfrow = c(2, 2))
plot.legend <- c("Weibull", "Gamma")
denscomp(list(fit.weibull, fit.gamma), fitcol = c("red", "blue"), legendtext = plot.legend)
qqcomp(list(fit.weibull, fit.gamma), fitcol = c("red", "blue"), legendtext = plot.legend)
cdfcomp(list(fit.weibull, fit.gamma), fitcol = c("red", "blue"), legendtext = plot.legend)
ppcomp(list(fit.weibull, fit.gamma), fitcol = c("red", "blue"), legendtext = plot.legend)

@NickCox indique à juste titre que le QQ-Plot (panneau supérieur droit) est le meilleur graphique unique pour juger et comparer les ajustements. Les densités ajustées sont difficiles à comparer. J'inclus également les autres graphiques par souci d'exhaustivité.

# Compare goodness of fit

gofstat(list(fit.weibull, fit.gamma))

Goodness-of-fit statistics
                             1-mle-weibull 2-mle-gamma
Kolmogorov-Smirnov statistic    0.06863193   0.1204876
Cramer-von Mises statistic      0.05673634   0.2060789
Anderson-Darling statistic      0.38619340   1.2031051

Goodness-of-fit criteria
                               1-mle-weibull 2-mle-gamma
Aikake's Information Criterion      519.8537    531.5180
Bayesian Information Criterion      524.5151    536.1795