Bagaimana cara menggambar grafik yang pas dan grafik sebenarnya dari distribusi gamma dalam satu plot?

Muat paket yang dibutuhkan.

library(ggplot2)
library(MASS)

Hasilkan 10.000 angka yang pas untuk distribusi gamma.

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

Gambarkan fungsi kerapatan probabilitas, seandainya kita tidak tahu distribusi x mana yang cocok.

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

Dari grafik, kita dapat belajar bahwa distribusi x cukup seperti distribusi gamma, jadi kami gunakan fitdistr()dalam paket MASSuntuk mendapatkan parameter bentuk dan laju distribusi gamma.

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

Gambar titik aktual (titik hitam) dan grafik pas (garis merah) di plot yang sama, dan inilah pertanyaannya, silakan lihat plotnya terlebih dahulu.

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

Saya punya dua pertanyaan:

1. Parameter nyata shape=2, rate=0.2dan parameter saya menggunakan fungsi fitdistr()untuk mendapatkan yang shape=2.01, rate=0.20. Keduanya hampir sama, tetapi mengapa grafik pas tidak cocok dengan titik sebenarnya dengan baik, pasti ada sesuatu yang salah dalam grafik pas, atau cara saya menggambar grafik pas dan poin aktual benar-benar salah, apa yang harus saya lakukan ?

2. Setelah saya mendapatkan parameter dari model yang saya buat, dengan cara apa saya mengevaluasi model, sesuatu seperti RSS (residual square sum) untuk model linier, atau nilai p shapiro.test(), ks.test()dan tes lainnya?

Saya miskin dalam pengetahuan statistik, bisakah Anda membantu saya?

ps: Saya sering mencari di Google, stackoverflow dan CV, tetapi tidak menemukan apa pun yang terkait dengan masalah ini

pertanyaan 1

Cara Anda menghitung kepadatan dengan tangan tampaknya salah. Tidak perlu membulatkan angka acak dari distribusi gamma. Seperti yang dicatat @Pascal, Anda dapat menggunakan histogram untuk memplot kerapatan poin. Dalam contoh di bawah ini, saya menggunakan fungsi densityuntuk memperkirakan kerapatan dan plot sebagai titik. Saya menyajikan kecocokan dengan poin dan histogram:

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

Inilah solusi yang disediakan @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()

Pertanyaan 2

Untuk menilai kebaikan yang cocok saya sarankan paket fitdistrplus. Berikut adalah bagaimana ini dapat digunakan untuk mencocokkan dua distribusi dan membandingkan kesesuaiannya secara grafis dan numerik. Perintah ini gofstatmencetak beberapa langkah, seperti AIC, BIC dan beberapa statistik gof seperti KS-Test dll. Ini terutama digunakan untuk membandingkan kesesuaian distribusi yang berbeda (dalam hal ini gamma versus Weibull). Informasi lebih lanjut dapat ditemukan dalam jawaban saya di sini :

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 dengan tepat menyarankan bahwa QQ-Plot (panel kanan atas) adalah grafik tunggal terbaik untuk menilai dan membandingkan kecocokan. Kepadatan yang dipasang sulit untuk dibandingkan. Saya menyertakan grafis lainnya juga demi kelengkapan.

# 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