Wie zeichnet man einen angepassten Graphen und einen tatsächlichen Graphen der Gammaverteilung in einem Diagramm?

10

Laden Sie das benötigte Paket.

library(ggplot2)
library(MASS)

Generieren Sie 10.000 Zahlen, die an die Gammaverteilung angepasst sind.

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

Zeichnen Sie die Wahrscheinlichkeitsdichtefunktion, vorausgesetzt, wir wissen nicht, an welche Verteilung x angepasst ist.

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

pdf

Aus dem Diagramm können wir lernen, dass die Verteilung von x der Gammaverteilung sehr ähnlich ist. Daher verwenden wir sie fitdistr()im Paket MASS, um die Parameter für Form und Geschwindigkeit der Gammaverteilung zu erhalten.

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

Zeichnen Sie den tatsächlichen Punkt (schwarzer Punkt) und das angepasste Diagramm (rote Linie) im selben Diagramm. Hier ist die Frage: Schauen Sie sich zuerst das Diagramm an.

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

angepasste Grafik

Ich habe zwei Fragen:

  1. Die realen Parameter sind shape=2, rate=0.2und die Parameter, mit denen ich die Funktion fitdistr()abrufe, sind shape=2.01, rate=0.20. Diese beiden sind fast gleich, aber warum das angepasste Diagramm nicht gut zum tatsächlichen Punkt passt, muss etwas im angepassten Diagramm nicht stimmen, oder die Art und Weise, wie ich das angepasste Diagramm und die tatsächlichen Punkte zeichne, ist völlig falsch. Was soll ich tun? ?

  2. Nachdem ich die Parameter des Modells erhalte ich schaffen, in welcher Weise ich das Modell zu bewerten, so etwas wie RSS (Restquadratsumme) für lineares Modell oder den p-Wert shapiro.test(), ks.test()und anderen Test?

Ich bin arm an statistischen Kenntnissen. Könnten Sie mir bitte helfen?

ps: Ich habe viele Male in Google, Stackoverflow und CV gesucht, aber nichts im Zusammenhang mit diesem Problem gefunden

Ling Zhang
quelle
1
Ich habe diese Frage zuerst im Stackoverflow gestellt, aber es schien, dass diese Frage zum Lebenslauf gehört. Der Freund sagte, ich habe die Wahrscheinlichkeitsmassenfunktion und die Wahrscheinlichkeitsdichtefunktion falsch verstanden. Ich konnte sie nicht vollständig erfassen. Verzeihen Sie mir, dass ich diese Frage in erneut beantwortet habe CV
Ling Zhang
1
Ihre Dichteberechnung ist falsch. Eine einfache Methode zur Berechnung ist h <- hist(x, 1000, plot = FALSE); t1 <- data.frame(x = h$mids, y = h$density).
@Pascal du hast recht, ich habe Q1 gelöst, danke!
Ling Zhang
Siehe Antwort unten, densityFunktion ist nützlich.
Ich verstehe, nochmals vielen Dank für die Bearbeitung und Lösung meiner Frage
Ling Zhang

Antworten:

11

Frage 1

Die Art und Weise, wie Sie die Dichte von Hand berechnen, scheint falsch zu sein. Es ist nicht erforderlich, die Zufallszahlen aus der Gammaverteilung zu runden. Wie @Pascal feststellte, können Sie ein Histogramm verwenden, um die Dichte der Punkte zu zeichnen. Im folgenden Beispiel verwende ich die Funktion density, um die Dichte zu schätzen und als Punkte zu zeichnen. Ich präsentiere die Anpassung sowohl mit den Punkten als auch mit dem Histogramm:

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

Gammadichte

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

Anpassung der Gammadichte

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

Histogramm mit Passform

Hier ist die Lösung, die @Pascal bereitgestellt hat:

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

Histogrammdichtepunkte

Frage 2

Um die Passform zu beurteilen, empfehle ich das Paket fitdistrplus. Hier erfahren Sie, wie Sie zwei Verteilungen anpassen und ihre Anpassungen grafisch und numerisch vergleichen können. Der Befehl gofstatdruckt verschiedene Kennzahlen aus, z. B. AIC, BIC und einige Gof-Statistiken wie KS-Test usw. Diese werden hauptsächlich zum Vergleichen von Anpassungen verschiedener Verteilungen verwendet (in diesem Fall Gamma gegenüber Weibull). Weitere Informationen finden Sie in meiner Antwort hier :

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 weist zu Recht darauf hin, dass das QQ-Diagramm (oberes rechtes Feld) das beste einzelne Diagramm zum Beurteilen und Vergleichen von Anpassungen ist. Angepasste Dichten sind schwer zu vergleichen. Der Vollständigkeit halber füge ich auch die anderen Grafiken hinzu.

Passungen vergleichen

# 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
COOLSerdash
quelle
1
Ich kann nicht überarbeiten, aber Sie haben ein Problem mit dem Backtick für fitdistrplusund gofstatin Ihrem Ansewer
2
Einzeilige Empfehlung: Das Quantil-Quantil-Diagramm ist der beste Einzelgraph für diesen Zweck. Der Vergleich der beobachteten und angepassten Dichten ist schwierig. Beispielsweise ist es schwierig, systematische Abweichungen bei hohen Werten zu erkennen, die wissenschaftlich und praktisch oft sehr wichtig sind.
Nick Cox
1
Ich bin froh, dass wir uns einig sind. Das OP beginnt mit 10.000 Punkten. Viele Probleme beginnen mit weitaus weniger und dann kann es problematisch sein, eine gute Vorstellung von der Dichte zu bekommen.
Nick Cox
1
@LingZhang Um Anpassungen zu vergleichen, können Sie sich den Wert des AIC ansehen. Die Anpassung mit dem niedrigsten AIC wird bevorzugt. Ich bin auch nicht der Meinung, dass die Weibull- und Gamma-Verteilung im QQ-Plot ziemlich gleich sind. Die Punkte der Weibull-Anpassung liegen im Vergleich zur Gamma-Anpassung näher an der Linie, insbesondere an den Schwänzen. Entsprechend ist der AIC für die Weibull-Anpassung im Vergleich zur Gamma-Anpassung kleiner.
COOLSerdash
1
Gerader ist besser. Siehe auch stats.stackexchange.com/questions/111010/…. Die Prinzipien sind dieselben. Die systematische Abweichung von der Linearität ist ein Problem.
Nick Cox