#############################
## 1. Binomial data
#############################
## (empirical) data
# seed for reproducibility:
set.seed(20200306)
x <- rbinom(100, size=25, prob=.25)
## ML-estimate
MLEstimator(x, BinomFamily(size = 25))
#############################
## 2. Poisson data
#############################
## Example: Rutherford-Geiger (1910); cf. Feller~(1968), Section VI.7 (a)
x <- c(rep(0, 57), rep(1, 203), rep(2, 383), rep(3, 525), rep(4, 532),
rep(5, 408), rep(6, 273), rep(7, 139), rep(8, 45), rep(9, 27),
rep(10, 10), rep(11, 4), rep(12, 0), rep(13, 1), rep(14, 1))
## ML-estimate
MLEstimator(x, PoisFamily())
#############################
## 3. Normal (Gaussian) location and scale
#############################
## (empirical) data
# seed for reproducibility:
set.seed(20200306)
x <- rnorm(100)
## ML-estimate
MLEstimator(x, NormLocationScaleFamily())
## compare:
c(mean(x),sd(x))
#############################
## 4. Gamma model
#############################
## (empirical) data
# seed for reproducibility:
set.seed(20200306)
x <- rgamma(50, scale = 0.5, shape = 3)
## parametric family of probability measures
G <- GammaFamily(scale = 1, shape = 2)
## Maximum likelihood estimator
(res <- MLEstimator(x = x, ParamFamily = G))
## Asymptotic (CLT-based) confidence interval
confint(res)
## some profiling
par(mfrow=c(1,2))
plot(profile(res))
par(mfrow=c(1,1))
## implementation of ML-estimator of package MASS
require(MASS)
(res1 <- fitdistr(x, "gamma"))
## comparison
## shape
estimate(res)[2]
## rate
1/estimate(res)[1]
## minor differences due to the fact that by default, fitdistr uses
## BFGS, while we use Nelder-Mead instead
## log-likelihood
res1$loglik
## negative log-likelihood
criterion(res)
## explicitely transforming to
## MASS parametrization:
mtrafo <- function(x){
nms0 <- names(c(main(param(G)),nuisance(param(G))))
nms <- c("shape","rate")
fval0 <- c(x[2], 1/x[1])
names(fval0) <- nms
mat0 <- matrix( c(0, -1/x[1]^2, 1, 0), nrow = 2, ncol = 2,
dimnames = list(nms,nms0))
list(fval = fval0, mat = mat0)}
G2 <- G
trafo(G2) <- mtrafo
res2 <- MLEstimator(x = x, ParamFamily = G2)
old <- getdistrModOption("show.details")
distrModoptions("show.details" = "minimal")
res1
res2
## some profiling
par(mfrow=c(1,2))
plot(profile(res2))
par(mfrow=c(1,1))
#############################
## 5. Cauchy Location Scale model
#############################
(C <- CauchyLocationScaleFamily())
loc.true <- 1
scl.true <- 2
## (empirical) data
# seed for reproducibility:
set.seed(20200306)
x <- rcauchy(50, location = loc.true, scale = scl.true)
## Maximum likelihood estimator
(res <- MLEstimator(x = x, ParamFamily = C))
## Asymptotic (CLT-based) confidence interval
confint(res)
Run the code above in your browser using DataLab