Problem 46: Bowhead whale.

Author

Prof. Eric A. Suess

Problem 46: Bowhead whale.

library(readr)

whale <- read_csv("~/classes/2023-2024/Fall 2023/Stat640/website/_homework/mle/whale/whale.csv")

Rows: 210 Columns: 1
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (1): times

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(whale)

# A tibble: 6 × 1
  times
  <dbl>
1 0.225
2 0.741
3 0.282
4 1.15 
5 0.437
6 0.346

convert the data_frame to a vector

The fitdistrplus package expects a vector as input.

whale <- whale$times
whale

  [1] 0.22520 0.74070 0.28250 1.14940 0.43670 0.34600 1.72410 0.58820 0.34480
 [10] 0.39840 1.08700 2.04080 0.25510 0.36230 0.22730 0.42190 0.24150 1.19050
 [19] 0.38170 0.16670 0.22830 0.71940 0.35710 0.45870 1.01010 0.55560 0.50760
 [28] 0.22730 0.52080 0.45660 1.78570 1.02040 0.45050 0.30030 0.17760 0.33670
 [37] 4.00000 1.29870 0.20530 0.34970 0.75760 0.38020 0.16420 0.39680 2.85710
 [46] 0.54050 0.74630 0.35840 1.13158 0.75190 1.42860 1.14940 0.34360 3.70370
 [55] 0.22520 0.22030 0.21880 0.28900 1.26580 0.21460 0.30210 0.25250 0.45870
 [64] 0.23980 0.68030 1.36990 0.42370 0.40320 0.63690 0.45250 0.34480 0.31450
 [73] 0.25250 0.47170 0.32680 0.86960 0.22220 0.19120 0.19840 0.24210 0.36360
 [82] 0.32680 0.16750 0.47850 0.18800 0.40490 0.37310 0.33670 1.18182 0.35590
 [91] 2.22220 1.78570 0.20790 0.34970 0.30030 0.42920 0.17180 0.31060 0.54050
[100] 0.33110 0.23530 0.23200 0.21510 0.25380 0.34600 0.21230 0.21690 1.49250
[109] 2.27270 0.31650 0.30210 0.51810 1.51520 0.23420 0.64100 0.32890 0.36230
[118] 0.29590 0.14750 0.36360 0.25710 1.14940 0.29590 0.36900 0.18120 0.40650
[127] 0.13370 1.72410 0.59880 0.51280 0.40490 0.33110 0.21410 0.43860 0.20160
[136] 0.41320 0.29150 0.20160 0.23470 0.28990 0.24940 0.31060 0.17240 0.46300
[145] 2.50000 0.18280 0.56180 0.34010 2.63160 0.41670 0.35840 0.27860 0.30860
[154] 0.93460 0.34840 0.14510 0.21790 0.37310 1.13640 0.32680 0.20960 0.34480
[163] 0.21790 0.27780 0.35210 1.56250 0.30580 0.42740 0.47620 0.27100 3.22580
[172] 0.44640 0.66700 1.06380 0.94340 0.28650 0.30860 0.25380 0.74070 0.21190
[181] 0.34970 0.14660 0.21280 0.33440 0.16860 0.40650 4.34780 0.31850 0.32790
[190] 0.20280 0.35710 0.72460 0.58820 0.63690 0.31550 0.38020 2.56410 0.36360
[199] 0.30120 0.86960 1.49250 0.23640 0.24040 0.29070 0.59170 0.20200 0.28410
[208] 0.13370 0.61350 0.23700

library(fitdistrplus)

Loading required package: MASS

Loading required package: survival

plotdist(whale, histo = TRUE, demp = TRUE)

fit gamma using MME and MLE

fit.gamma.mme <- fitdist(whale, "gamma", method = "mme")
summary(fit.gamma.mme)

Fitting of the distribution ' gamma ' by matching moments 
Parameters : 
       estimate
shape 0.8002543
rate  1.3291210
Loglikelihood:  -116.5212   AIC:  237.0423   BIC:  243.7365

fit.gamma.mle <- fitdist(whale, "gamma", method = "mle")
summary(fit.gamma.mle)

Fitting of the distribution ' gamma ' by maximum likelihood 
Parameters : 
      estimate Std. Error
shape 1.611251  0.1439041
rate  2.676063  0.2797996
Loglikelihood:  -90.89708   AIC:  185.7942   BIC:  192.4884 
Correlation matrix:
          shape      rate
shape 1.0000000 0.8541993
rate  0.8541993 1.0000000

plot the log-likelihood

llplot(fit.gamma.mle)

llplot(fit.gamma.mle, expansion = 2)

compare

denscomp(list(fit.gamma.mme, fit.gamma.mle), legendtext = c("MME", "MLE"))

cdfcomp(list(fit.gamma.mme, fit.gamma.mle), legendtext = c("MME", "MLE"))

qqcomp(list(fit.gamma.mme, fit.gamma.mle), legendtext = c("MME", "MLE"))

bootstrap the MME

b1 <- bootdist(fit.gamma.mme, niter=101)
print(b1)

Parameter values obtained with parametric bootstrap 
      shape     rate
1 0.9418513 1.576284
2 0.9410039 1.483790
3 0.9691874 1.440645
4 0.9145934 1.487089
5 0.8713925 1.397855
6 0.9033110 1.539314

plot(b1)

hist(b1$estim$shape)

hist(b1$estim$rate)

summary(b1)

Parametric bootstrap medians and 95% percentile CI 
         Median      2.5%    97.5%
shape 0.8230687 0.6307503 1.023911
rate  1.3840386 1.0526418 1.907737

bootstrap the MLE

b2 <- bootdist(fit.gamma.mle, niter=101)
print(b2)

Parameter values obtained with parametric bootstrap 
     shape     rate
1 1.585549 2.596792
2 1.659788 2.648884
3 1.771854 3.195506
4 1.664727 2.713174
5 1.431631 2.490766
6 1.745494 2.760940

hist(b2$estim$shape)

hist(b2$estim$rate)

plot(b2)

summary(b2)

Parametric bootstrap medians and 95% percentile CI 
        Median     2.5%    97.5%
shape 1.628309 1.410624 1.905667
rate  2.722925 2.327627 3.254625