<- 100
n
<- .5; a.min <- 0.1; a.max <- 2;
alpha <- 2; l.max <- 10
lambda
# generate random data
<- rgamma(n, shape = alpha, rate = lambda)
x
<- x # save a copy of the original data x.orig
MLE Gamma
MME fitted gamma
\(X_1, X_2, ...,X_n\) iid \(Gamma(\alpha,\lambda)\), find the MMEs.
\(f(x|\alpha,\lambda)= \frac{\lambda^\alpha}{\Gamma(\alpha)} x^{\alpha-1} e^{-\lambda x}\)
make a histogram
use the density estimator since this is a continuous random variable
hist(x.orig, probability=T)
lines(density(x.orig))
MMEs
<- (1/n)*sum( (x - mean(x))^2 )
sigma2.hat
<- mean(x)^2/sigma2.hat
alpha.hat.mme alpha.hat.mme
[1] 0.3472268
<- mean(x)/sigma2.hat
lambda.hat.mme lambda.hat.mme
[1] 1.63452
The AIC and BIC are defined as follows, with \(n\) the sample size and \(k\) the number of parameters in the model.
\[AIC = -2\log(L(\hat{\theta})) + 2k\]
\[BIC = -2\log(L(\hat{\theta})) + k\log(n)\]
library(fitdistrplus)
Loading required package: MASS
Loading required package: survival
<- fitdist(x.orig, "gamma", method = "mme")
fit.gam summary(fit.gam)
Fitting of the distribution ' gamma ' by matching moments
Parameters :
estimate Std. Error
shape 0.3472268 0.9672571
rate 1.6345201 5.3316147
Loglikelihood: 76.36084 AIC: -148.7217 BIC: -143.5113
Correlation matrix:
shape rate
shape 1.0000000 0.8540047
rate 0.8540047 1.0000000
plot the fitted model to the simulated data
denscomp(list(fit.gam), legendtext = c("Gamma"))
cdfcomp(list(fit.gam), legendtext = c("Gamma"))
qqcomp(list(fit.gam), legendtext = c("Gamma"))
ppcomp(list(fit.gam), legendtext = c("Gamma"))
MLE fitted gamma
<- fitdist(x.orig, "gamma", method = "mle")
fit.gam summary(fit.gam)
Fitting of the distribution ' gamma ' by maximum likelihood
Parameters :
estimate Std. Error
shape 0.4788249 0.05567706
rate 2.2542167 0.41812618
Loglikelihood: 79.86981 AIC: -155.7396 BIC: -150.5293
Correlation matrix:
shape rate
shape 1.0000000 0.6268848
rate 0.6268848 1.0000000
plot the fitted model to the simulated data
denscomp(list(fit.gam), legendtext = c("Gamma"))
cdfcomp(list(fit.gam), legendtext = c("Gamma"))
qqcomp(list(fit.gam), legendtext = c("Gamma"))
ppcomp(list(fit.gam), legendtext = c("Gamma"))