Wine Project: Part B
Adam Deuber
10/1/2020
- Suppose the population mean of the variable “density” is µ , do the following inferences:
- Provide an estimate of µ based on the sample
library("readxl")
red_wine <- read_excel("C:/Users/Adam Deuber/OneDrive/UC/BANA Masters/Statistical Methods/Homework/Wine Project/winequality-red.xlsx")
attach(red_wine)
mean_density <- mean(density)
mean_density## [1] 0.9967467- Use the Central Limit Theorem (CLT) to quantify the variability of your estimate
standev_density <- sd(density)
standev_density## [1] 0.001887334hist(density)
As seen by these values, the means are very close to 1 and standard deviations are very close to 0. The Central Limit Theorem is showcased here as the density sample reflects a normal distribution. Additionally, the histogram further reflect the distribution for the population and sample.
- Use the CLT to give a 95% confidence interval for µ.
lower <- mean(density) - 2*(sd(density)/sqrt(length(density)))
upper <- mean(density) + 2*(sd(density)/sqrt(length(density)))
lower## [1] 0.9966523upper## [1] 0.9968411As seen above, the 95% interval falls between 0.9966523 to 0.9968411.
- Use the bootstrap method to do parts b and c, and compare the results with those obtained from the CLT. State your findings.
density_mu.set <- NULL
for(k in 1:3500) {
bootstrap <- sample(density, size=1599, replace=T)
mu_density <- mean(bootstrap)
density_mu.set[k] <- mu_density
}
mean(density_mu.set)## [1] 0.9967449sd(density_mu.set)## [1] 4.687592e-05conf_quantile <- quantile(density_mu.set, probs = c(0.025,0.975))
conf_quantile## 2.5% 97.5%
## 0.9966546 0.9968377Utilizing the bootstrap method, we can see that the results come to very similar answers as the CLT. This for loop proves both ways are effective to prove the normal distribution of density. The sample will provide very slight differences eacht time the for loop is ran, but the results are always close to the CLT.
- Can we use a normal distribution to model “density”? If yes, what are the maximum likelihood estimates of the mean and standard deviation? Please provide their standard errors as well.
library(stats4)
hist(density)
density_10 <- density*10
density_10## [1] 9.9780 9.9680 9.9700 9.9800 9.9780 9.9780 9.9640 9.9460 9.9680
## [10] 9.9780 9.9590 9.9780 9.9430 9.9740 9.9860 9.9860 9.9690 9.9680
## [19] 9.9740 9.9690 9.9680 9.9820 9.9660 9.9680 9.9680 9.9550 9.9620
## [28] 9.9660 9.9720 9.9640 9.9580 9.9660 9.9660 9.9930 9.9570 9.9860
## [37] 9.9750 9.9680 9.9400 9.9780 9.9780 9.9760 9.9680 9.9680 9.9620
## [46] 9.9340 9.9700 9.9690 9.9580 9.9540 9.9710 9.9560 9.9550 9.9700
## [55] 9.9550 9.9780 9.9710 9.9830 9.9750 9.9620 9.9800 9.9680 9.9680
## [64] 9.9660 9.9620 9.9620 9.9680 9.9620 9.9690 9.9620 9.9670 9.9620
## [73] 9.9610 9.9760 9.9840 9.9860 9.9860 9.9660 9.9580 9.9720 9.9580
## [82] 9.9740 9.9700 9.9690 9.9590 9.9610 9.9720 9.9660 9.9780 9.9780
## [91] 9.9640 9.9720 9.9720 9.9660 9.9380 9.9320 9.9650 9.9630 9.9670
## [100] 9.9720 9.9720 9.9590 9.9720 9.9670 9.9600 9.9670 9.9730 9.9690
## [109] 9.9760 9.9690 9.9640 9.9640 9.9640 9.9880 9.9640 9.9660 9.9780
## [118] 9.9700 9.9620 9.9590 9.9620 9.9620 9.9800 9.9760 9.9600 9.9840
## [127] 9.9640 9.9640 9.9620 9.9550 9.9680 9.9370 9.9370 9.9520 9.9590
## [136] 9.9650 9.9680 9.9720 9.9620 9.9620 9.9650 9.9680 9.9160 9.9540
## [145] 9.9160 9.9680 9.9440 9.9680 9.9590 9.9750 9.9580 9.9960 9.9500
## [154] 9.9500 9.9730 9.9730 9.9730 9.9730 9.9690 9.9680 9.9590 9.9650
## [163] 9.9640 9.9820 9.9810 9.9610 9.9560 9.9560 9.9530 9.9640 9.9720
## [172] 9.9720 9.9720 9.9610 9.9610 9.9580 9.9610 9.9720 9.9690 9.9690
## [181] 9.9690 9.9720 9.9720 9.9710 9.9700 9.9710 9.9720 9.9760 9.9680
## [190] 9.9680 9.9760 9.9740 9.9690 9.9640 9.9640 9.9680 9.9680 9.9810
## [199] 9.9240 9.9480 9.9695 9.9750 9.9545 9.9610 9.9615 9.9940 9.9940
## [208] 9.9625 9.9660 9.9800 9.9585 9.9685 9.9800 9.9590 9.9740 9.9740
## [217] 9.9620 9.9670 9.9570 9.9655 9.9640 9.9570 9.9525 9.9815 9.9745
## [226] 9.9615 9.9640 9.9840 9.9615 9.9685 9.9270 9.9625 9.9800 9.9685
## [235] 9.9685 9.9675 9.9675 9.9675 9.9675 9.9685 9.9710 9.9925 9.9565
## [244] 10.0005 10.0005 9.9830 9.9685 9.9630 9.9615 9.9830 9.9850 9.9720
## [253] 9.9860 9.9590 9.9720 9.9675 9.9800 9.9680 9.9680 9.9965 9.9625
## [262] 9.9565 9.9575 9.9630 9.9990 9.9680 10.0025 9.9730 9.9870 9.9960
## [271] 9.9650 9.9960 9.9935 9.9735 9.9915 9.9650 9.9870 9.9960 9.9760
## [280] 9.9910 9.9860 9.9720 9.9710 9.9910 10.0015 10.0015 9.9970 9.9670
## [289] 9.9760 10.0100 9.9760 9.9790 9.9940 9.9640 10.0140 10.0010 9.9810
## [298] 9.9855 9.9930 9.9845 9.9650 9.9730 9.9940 9.9600 9.9730 9.9980
## [307] 9.9660 9.9940 9.9940 9.9680 9.9980 9.9600 9.9815 9.9760 9.9645
## [316] 9.9620 9.9960 9.9865 9.9890 9.9865 9.9890 9.9975 9.9735 9.9900
## [325] 10.0150 10.0150 9.9780 9.9800 10.0020 9.9810 9.9800 9.9800 9.9730
## [334] 9.9580 9.9630 9.9920 9.9480 9.9740 10.0080 9.9970 10.0000 9.9720
## [343] 9.9820 9.9820 10.0060 9.9790 9.9550 9.9860 9.9790 9.9940 10.0040
## [352] 9.9940 9.9740 10.0180 9.9120 9.9550 10.0010 9.9700 10.0000 10.0100
## [361] 9.9730 9.9820 9.9920 9.9880 10.0220 9.9790 10.0220 10.0000 9.9980
## [370] 9.9620 9.9750 9.9670 9.9670 9.9590 10.0140 9.9840 10.0030 9.9620
## [379] 9.9880 9.9790 9.9740 10.0140 9.9740 9.9740 9.9800 9.9790 9.9810
## [388] 9.9720 9.9810 9.9820 9.9240 10.0140 9.9800 9.9700 9.9970 10.0040
## [397] 9.9940 9.9880 9.9880 9.9630 9.9940 9.9490 10.0000 9.9940 9.9740
## [406] 9.9660 9.9810 9.9960 9.9700 9.9810 9.9870 9.9870 9.9660 9.9860
## [415] 9.9820 10.0140 9.9700 9.9560 9.9960 9.9760 9.9720 9.9520 9.9550
## [424] 9.9780 9.9550 9.9520 9.9490 9.9880 9.9780 9.9970 9.9780 9.9740
## [433] 9.9580 10.0040 9.9960 10.0040 9.9580 9.9860 9.9960 9.9510 9.9870
## [442] 9.9880 10.0320 9.9630 9.9340 9.9860 9.9960 9.9680 9.9670 9.9940
## [451] 9.9940 9.9790 9.9680 9.9980 9.9520 9.9880 9.9860 9.9880 9.9980
## [460] 10.0080 9.9760 9.9820 9.9700 9.9720 9.9870 10.0060 9.9820 9.9470
## [469] 10.0000 9.9520 9.9960 9.9700 9.9950 9.9710 9.9820 9.9880 9.9800
## [478] 9.9770 9.9880 9.9840 10.0260 9.9640 9.9860 9.9860 9.9700 9.9760
## [487] 9.9760 9.9820 10.0020 9.9840 9.9840 9.9520 9.9500 10.0020 9.9520
## [496] 9.9720 9.9840 9.9660 9.9720 10.0020 9.9840 9.9900 9.9900 9.9760
## [505] 9.9760 9.9680 9.9760 10.0000 9.9940 9.9860 9.9780 9.9940 9.9800
## [514] 9.9730 9.9730 10.0100 10.0100 9.9940 9.9720 9.9660 9.9820 9.9800
## [523] 9.9790 9.9780 9.9760 9.9980 9.9660 9.9620 9.9880 9.9740 9.9900
## [532] 10.0040 10.0040 9.9560 10.0000 9.9900 9.9740 9.9720 10.0140 9.9960
## [541] 9.9880 9.9880 9.9660 9.9930 10.0080 9.9820 9.9830 9.9870 9.9940
## [550] 9.9750 9.9580 9.9800 9.9800 9.9340 10.0315 10.0315 10.0020 10.0315
## [559] 10.0020 10.0210 9.9940 9.9780 9.9820 9.9710 10.0210 9.9940 9.9910
## [568] 9.9910 9.9900 9.9460 10.0030 9.9460 9.9780 9.9910 9.9920 9.9930
## [577] 9.9780 9.9820 9.9820 9.9760 10.0020 10.0020 10.0000 9.9760 10.0020
## [586] 9.9800 9.9860 9.9670 9.9170 9.9760 9.9870 9.9220 9.9870 10.0040
## [595] 9.9780 9.9980 9.9940 9.9800 9.9670 10.0000 9.9620 10.0060 9.9750
## [604] 10.0060 9.9840 9.9840 9.9730 9.9820 10.0260 9.9560 9.9830 10.0060
## [613] 9.9650 9.9620 9.9690 9.9970 9.9970 9.9820 10.0060 9.9880 9.9760
## [622] 9.9760 9.9740 9.9740 9.9970 9.9970 9.9880 9.9880 9.9760 9.9640
## [631] 9.9760 9.9880 9.9640 10.0100 9.9640 9.9540 9.9790 9.9780 9.9640
## [640] 9.9700 9.9910 9.9920 9.9910 9.9920 9.9910 9.9840 9.9720 9.9670
## [649] 9.9640 9.9480 9.9860 10.0100 9.9760 9.9760 9.9790 10.0040 9.9860
## [658] 9.9800 9.9730 9.9700 9.9730 9.9780 9.9720 9.9740 10.0000 9.9810
## [667] 9.9800 9.9880 9.9660 9.9880 9.9590 9.9680 9.9800 9.9680 9.9840
## [676] 9.9840 9.9840 9.9790 9.9830 9.9720 10.0040 9.9620 9.9800 9.9700
## [685] 9.9800 9.9700 9.9700 9.9980 9.9620 9.9760 9.9660 9.9980 9.9790
## [694] 9.9840 9.9840 9.9210 9.9720 9.9720 10.0010 9.9860 9.9840 9.9720
## [703] 9.9700 9.9720 9.9800 9.9900 9.9560 9.9660 9.9600 9.9830 9.9850
## [712] 9.9980 9.9920 9.9760 9.9800 9.9670 9.9760 9.9680 9.9760 9.9760
## [721] 9.9760 9.9820 9.9680 9.9650 9.9580 9.9760 9.9940 9.9700 9.9700
## [730] 9.9500 9.9820 9.9940 9.9660 9.9800 9.9870 9.9630 9.9630 10.0020
## [739] 9.9830 9.9900 9.9650 9.9788 9.9720 10.0024 10.0010 9.9768 9.9780
## [748] 9.9782 9.9761 9.9768 9.9803 9.9803 9.9785 9.9803 9.9656 9.9525
## [757] 9.9488 9.9656 9.9656 9.9823 9.9779 9.9738 9.9701 9.9738 9.9888
## [766] 9.9888 9.9738 9.9938 9.9744 9.9668 9.9744 9.9780 9.9782 9.9730
## [775] 9.9727 9.9586 9.9612 9.9610 9.9788 9.9745 9.9676 9.9732 9.9814
## [784] 9.9732 9.9746 9.9820 9.9820 9.9910 9.9910 9.9800 9.9708 9.9818
## [793] 9.9745 9.9639 9.9531 9.9786 9.9746 9.9526 9.9870 9.9870 9.9641
## [802] 9.9735 9.9264 9.9710 9.9682 9.9356 9.9386 9.9356 9.9702 9.9693
## [811] 9.9562 10.0012 9.9818 9.9462 9.9939 9.9818 9.9840 9.9545 9.9632
## [820] 9.9976 9.9606 9.9154 9.9730 9.9730 9.9682 9.9624 9.9510 9.9624
## [829] 9.9417 9.9376 9.9632 9.9376 9.9832 9.9836 9.9694 9.9655 9.9064
## [838] 9.9064 9.9672 9.9647 9.9736 9.9629 9.9708 9.9630 9.9590 9.9689
## [847] 9.9689 9.9770 9.9689 9.9689 9.9708 9.9708 9.9801 9.9652 9.9652
## [856] 9.9538 9.9652 9.9594 9.9686 9.9438 9.9746 9.9357 9.9628 9.9746
## [865] 9.9746 9.9748 9.9438 9.9438 9.9438 9.9578 9.9371 9.9522 9.9576
## [874] 9.9552 9.9664 9.9614 9.9517 9.9371 9.9787 9.9745 9.9576 9.9533
## [883] 9.9536 9.9745 9.9787 9.9824 9.9836 9.9577 9.9491 10.0289 9.9576
## [892] 9.9743 9.9774 9.9743 9.9745 9.9550 9.9444 9.9550 9.9444 9.9892
## [901] 9.9562 9.9736 9.9736 9.9620 9.9620 9.9800 9.9640 9.9480 9.9528
## [910] 9.9331 9.9577 9.9901 9.9674 9.9639 9.9331 9.9512 9.9395 9.9824
## [919] 9.9640 9.9504 9.9786 9.9640 9.9504 9.9824 9.9516 9.9604 9.9786
## [928] 9.9736 9.9516 9.9468 9.9746 9.9748 9.9710 9.9748 9.9746 9.9543
## [937] 9.9543 9.9791 9.9356 9.9425 9.9509 9.9484 9.9834 9.9864 9.9498
## [946] 9.9566 9.9745 9.9408 9.9552 9.9552 9.9552 9.9408 9.9536 9.9458
## [955] 9.9538 9.9648 9.9568 9.9613 9.9519 9.9735 9.9518 9.9592 9.9654
## [964] 9.9546 9.9518 9.9554 9.9604 9.9733 9.9430 9.9669 9.9724 9.9724
## [973] 9.9643 9.9605 9.9658 9.9700 9.9700 9.9801 9.9416 9.9690 9.9774
## [982] 9.9712 9.9418 9.9774 9.9690 9.9562 9.9470 9.9693 9.9596 9.9556
## [991] 9.9596 9.9694 9.9554 9.9694 9.9918 9.9697 9.9378 9.9378 9.9554
## [1000] 9.9162 9.9495 9.9676 9.9516 9.9280 9.9603 9.9280 9.9516 9.9516
## [1009] 9.9549 9.9722 9.9354 9.9570 9.9604 9.9635 9.9454 9.9598 9.9486
## [1018] 9.9007 9.9007 9.9636 9.9642 9.9642 9.9584 9.9506 9.9568 9.9822
## [1027] 9.9364 9.9378 9.9586 9.9568 9.9488 9.9514 9.9854 9.9592 9.9739
## [1036] 9.9683 9.9356 9.9672 9.9530 9.9692 9.9756 9.9547 9.9692 9.9859
## [1045] 9.9294 9.9438 9.9612 9.9634 9.9702 9.9704 9.9634 9.9702 9.9258
## [1054] 9.9426 9.9747 9.9747 9.9586 9.9784 9.9710 9.9586 9.9810 9.9462
## [1063] 9.9560 9.9565 9.9418 9.9630 9.9358 9.9572 9.9572 9.9700 9.9498
## [1072] 9.9892 9.9680 9.9720 9.9892 9.9769 9.9534 9.9817 9.9817 9.9316
## [1081] 9.9471 9.9316 9.9685 9.9617 9.9685 9.9529 9.9738 9.9400 9.9735
## [1090] 9.9735 9.9451 9.9546 9.9479 9.9615 9.9655 9.9772 9.9655 9.9666
## [1099] 9.9358 9.9666 9.9392 9.9388 9.9402 9.9388 9.9360 9.9374 9.9376
## [1108] 9.9523 9.9855 9.9820 9.9593 9.9471 9.9396 9.9698 9.9020 9.9572
## [1117] 9.9572 9.9572 9.9252 9.9256 9.9235 9.9352 9.9220 9.9674 9.9557
## [1126] 9.9394 9.9150 9.9379 9.9840 9.9798 9.9770 9.9341 9.9330 9.9684
## [1135] 9.9524 9.9460 9.9774 9.9774 9.9786 9.9764 9.9690 9.9624 9.9354
## [1144] 9.9672 9.9588 9.9672 9.9600 9.9702 9.9526 9.9585 9.9473 9.9352
## [1153] 9.9616 9.9622 9.9544 9.9616 9.9524 9.9240 9.9572 9.9728 9.9676
## [1162] 9.9551 9.9434 9.9692 9.9692 9.9634 9.9709 9.9528 9.9384 9.9502
## [1171] 9.9667 9.9522 9.9649 9.9716 9.9716 9.9541 9.9572 9.9318 9.9346
## [1180] 9.9599 9.9599 9.9478 9.9754 9.9572 9.9538 9.9346 9.9630 9.9346
## [1189] 9.9538 9.9616 9.9652 9.9551 9.9439 9.9551 9.9588 9.9633 9.9686
## [1198] 9.9458 9.9510 9.9686 9.9458 9.9419 9.9516 9.9878 9.9534 9.9534
## [1207] 9.9534 9.9752 9.9534 9.9428 9.9498 9.9570 9.9498 9.9659 9.9628
## [1216] 9.9488 9.9677 9.9408 9.9600 9.9478 9.9734 9.9734 9.9678 9.9638
## [1225] 9.9920 9.9922 9.9592 9.9769 9.9157 9.9718 9.9470 9.9621 9.9718
## [1234] 9.9614 9.9242 9.9572 9.9659 9.9242 9.9798 9.9488 9.9494 9.9729
## [1243] 9.9414 9.9830 9.9721 9.9655 9.9502 9.9655 9.9490 9.9514 9.9514
## [1252] 9.9630 9.9627 9.9569 9.9628 9.9358 9.9499 9.9633 9.9538 9.9538
## [1261] 9.9632 9.9437 9.9726 9.9780 9.9456 9.9564 9.9564 9.9600 9.9668
## [1270] 9.9080 9.9084 9.9350 9.9385 9.9494 9.9440 9.9688 9.9566 9.9636
## [1279] 9.9688 9.9480 9.9560 9.9560 9.9619 9.9734 9.9476 9.9734 9.9328
## [1288] 9.9286 9.9914 9.9914 9.9521 9.9638 9.9362 9.9672 9.9638 9.9558
## [1297] 9.9558 9.9323 9.9191 9.9476 9.9444 9.9598 9.9576 9.9501 9.9550
## [1306] 9.9612 9.9790 9.9524 9.9790 9.9600 9.9612 9.9290 9.9532 9.9558
## [1315] 9.9556 9.9616 9.9258 9.9638 9.9616 9.9652 9.9796 9.9258 9.9392
## [1324] 9.9495 9.9480 9.9480 9.9480 9.9480 9.9550 9.9616 9.9616 9.9668
## [1333] 9.9581 9.9760 9.9608 9.9387 9.9448 9.9448 9.9448 9.9538 9.9538
## [1342] 9.9538 9.9589 9.9538 9.9852 9.9613 9.9472 9.9587 9.9587 9.9518
## [1351] 9.9654 9.9332 9.9690 9.9690 9.9568 9.9464 9.9464 9.9648 9.9736
## [1360] 9.9704 9.9699 9.9724 9.9704 9.9652 9.9500 9.9710 9.9576 9.9725
## [1369] 9.9656 9.9396 9.9623 9.9476 9.9623 9.9664 9.9471 9.9586 9.9460
## [1378] 9.9530 9.9609 9.9557 9.9557 9.9534 9.9613 9.9613 9.9593 9.9668
## [1387] 9.9610 9.9610 9.9718 9.9524 9.9292 9.9612 9.9420 9.9612 9.9669
## [1396] 9.9745 9.9745 9.9584 9.9660 9.9560 9.9620 9.9620 9.9630 9.9600
## [1405] 9.9850 9.9520 9.9740 9.9500 9.9500 9.9500 9.9360 9.9630 9.9740
## [1414] 9.9630 9.9940 9.9660 9.9940 9.9520 9.9500 9.9560 9.9500 9.9650
## [1423] 9.9560 9.9680 9.9740 9.9740 9.9710 9.9570 9.9620 9.9600 9.9560
## [1432] 9.9580 9.9440 9.9550 10.0369 10.0369 9.9914 9.9577 9.9600 9.9566
## [1441] 9.9470 9.9713 9.9712 9.9322 9.9566 9.9713 9.9712 9.9701 9.9683
## [1450] 9.9472 9.9470 9.9732 9.9374 9.9706 9.9974 9.9467 9.9236 9.9706
## [1459] 9.9603 9.9458 9.9724 9.9664 9.9545 9.9522 9.9580 9.9580 9.9728
## [1468] 9.9648 9.9728 9.9705 9.9578 9.9334 9.9656 9.9336 10.0242 9.9182
## [1477] 10.0242 9.9182 9.9808 9.9828 9.9498 9.9828 9.9719 9.9542 9.9592
## [1486] 9.9606 9.9546 9.9496 9.9420 9.9448 9.9344 9.9420 9.9348 9.9636
## [1495] 9.9459 9.9492 9.9636 9.9508 9.9582 9.9508 9.9508 9.9642 9.9638
## [1504] 9.9555 9.9605 9.9600 9.9562 9.9605 9.9814 9.9410 9.9661 9.9712
## [1513] 9.9588 9.9294 9.9842 9.9842 9.9633 9.9489 9.9647 9.9665 9.9489
## [1522] 9.9585 9.9633 9.9530 9.9522 9.9552 9.9553 9.9510 9.9714 9.9608
## [1531] 9.9444 9.9631 9.9573 9.9717 9.9397 9.9590 9.9528 9.9484 9.9538
## [1540] 9.9714 9.9646 9.9666 9.9472 9.9758 9.9550 9.9531 9.9575 9.9306
## [1549] 9.9783 9.9419 9.9783 9.9768 9.9586 9.9765 9.9627 9.9514 9.9636
## [1558] 9.9627 9.9787 9.9622 9.9622 9.9622 9.9546 9.9546 9.9546 9.9489
## [1567] 9.9494 9.9546 9.9629 9.9396 9.9340 9.9514 9.9632 9.9467 9.9677
## [1576] 9.9474 9.9588 9.9622 9.9540 9.9402 9.9470 9.9402 9.9362 9.9578
## [1585] 9.9484 9.9494 9.9492 9.9483 9.9414 9.9770 9.9314 9.9402 9.9574
## [1594] 9.9651 9.9490 9.9512 9.9574 9.9547 9.9549sd_density_10 <- sd(density_10)
dnorm(density_10[1], mean=mean(density_10), sd=sd(density_10))## [1] 18.08945minuslog.lik<-function(mu, sigma){
log.lik<-0
for(i in 1:1599) {
log.lik<-log.lik+log(dnorm(density_10[i], mean=mu,
sd=sigma))
}
return(-log.lik)
}
minuslog.lik(80, 20)## [1] 16062.54minuslog.lik(100, 30)## [1] 14108.6minuslog.lik(0, 1)## [1] 80900.31est <- stats4::mle(minuslog = minuslog.lik,
start = list(mu = mean(density_10),
sigma = sd(density_10)),
lower=c(0, 0))
summary(est)## Maximum likelihood estimation
##
## Call:
## stats4::mle(minuslogl = minuslog.lik, start = list(mu = mean(density_10),
## sigma = sd(density_10)), lower = c(0, 0))
##
## Coefficients:
## Estimate Std. Error
## mu 9.96746679 0.0004719810
## sigma 0.01887334 0.0003296881
##
## -2 log L: -8159.31Yes, we can use a normal distribution. The maximum likelihood estimates are 9.976 for the mean and 0.018 for the standard deviation. The standard error for the mean is 0.000472 for the mean and 0.000330 for the standard deviation.
- Suppose the population mean of the variable “residual sugar” is µ , answer the following questions.
- Provide an estimate of µ based on the sample
mean_resid_sugar <- mean(`residual sugar`)
mean_resid_sugar## [1] 2.538806- Noting that the sample distribution of “residual sugar” is highly skewed, can we use the CLT to quantify the variability of your estimate? Can we use the CLT to give a 95% confidence interval for µ? If yes, please give your solution. If no, explain why.
standev_sugar <- sd(`residual sugar`)
standev_sugar## [1] 1.409928hist(`residual sugar`)
lower_sug <- mean(`residual sugar`) - 2*(sd(`residual sugar`)/sqrt(length(`residual sugar`)))
upper_sug <- mean(`residual sugar`) + 2*(sd(`residual sugar`)/sqrt(length(`residual sugar`)))
lower_sug## [1] 2.468287upper_sug## [1] 2.609324Yes, we can use the CLT to quantify the variability of residual sugar. In addition we can use it to give a 95% confidence interval. As shown in the code, the bounds come to 2.468 and 2.609.
- Use the bootstrap method to do part b. Is the bootstrap confidence interval symmetric?
mu_sugar.set <- NULL
for(k in 1:3500) {
sugar.bootstrap <- sample(`residual sugar`, size=1599, replace=T)
mu_sugar <- mean(sugar.bootstrap)
mu_sugar.set[k] <- mu_sugar
}
mean(mu_sugar.set)## [1] 2.538663sd(mu_sugar.set)## [1] 0.03615705conf_quantile <- quantile(mu_sugar.set, probs = c(0.025,0.975))
conf_quantile## 2.5% 97.5%
## 2.46977 2.61076hist(mu_sugar.set, freq = FALSE)
lines(density(mu_sugar.set), lwd=5, col='blue')
Yes, we can use the bootstrap method as well. As seen by the histrogram, the distribution is mostly symmetric.
- Can we use a normal distribution to model “density”? If no, what distribution do you think can approximate its empirical distribution? What parameters are needed to characterize such a distribution? what are their maximum likelihood estimates? Please provide their standard errors as well.
minuslog.lik <- function(mu, sigma) {
log.lik <- 0
for(i in 1:1599) {
log.lik = log.lik + log(dlnorm(x=`residual sugar`[i], meanlog=mu, sdlog=sigma))
}
return(-log.lik)
}
estimate_sugar <- mle(minuslog=minuslog.lik,
start = list(mu=log(mean(`residual sugar`)), sigma=log(sd(`residual sugar`))),lower=c(0,0))
summary(estimate_sugar)## Maximum likelihood estimation
##
## Call:
## mle(minuslogl = minuslog.lik, start = list(mu = log(mean(`residual sugar`)),
## sigma = log(sd(`residual sugar`))), lower = c(0, 0))
##
## Coefficients:
## Estimate Std. Error
## mu 0.8502316 0.008936108
## sigma 0.3573326 0.006318587
##
## -2 log L: 3965.773No, we cannot use a normal distribution to model residual sugar. The histogram is skewed to the right as seen above. As such, I used lognormal fucntions to express the distribution. Mean and standard deviation are needed to characterize this distribution.
The maximum likelihood estimates are 0.850 for the mean and 0.357 for the standard deviation. The standard error for the mean is 0.008936 for the mean and 0.006318 for the standard deviation.
- We classify those wines as “excellent” if their rating is at least 7. Suppose the population proportion of excellent wines is p. Do the following: proportion of excellent wines is p. Do the following:
- Use the CLT to derive a 95% confidence interval for p;
red_wine$excellent <- as.numeric(red_wine$quality > 6)
hat <- mean(red_wine$excellent)
variation_excel <- sqrt(hat*(1 - hat) / length(red_wine$excellent))
variation_excel## [1] 0.008564681lower_qual <- hat - 2*(variation_excel)
upper_qual <- hat + 2*(variation_excel)
lower_qual## [1] 0.1185805upper_qual## [1] 0.1528392The 95% confidence interval gives us bounds of 0.1186 and 0.1528.
- Use the bootstrap method to derive a 95% confidence interval for p;
mu_quality.set <- NULL
for(k in 1:3500) {
quality.bootstrap <- sample(red_wine$excellent, size=1599, replace=T)
mu_quality <- mean(quality.bootstrap)
mu_quality.set[k] <- mu_quality
}
mean(mu_quality.set)## [1] 0.1358249sd(mu_quality.set)## [1] 0.008656594conf_quantile_qual <- quantile(mu_quality.set, probs = c(0.025,0.975))
conf_quantile_qual## 2.5% 97.5%
## 0.1194497 0.1532208Bootstrap gives us a confidence interval that is similar to the interval seen above for every time the loop is run.
- Compare the two intervals. Is there any difference worth our attention?
The intervals are mostly the same and do not require further attention.
- What is the maximum likelihood estimate of p and its standard error?
minuslog.lik<-function(p){
log.lik<-0
for(i in 1:1599) {
log.lik<-log.lik+log(dbinom(red_wine$excellent[i], size = 1, prob = p))
}
return(-log.lik)
}
est <- stats4::mle(minuslog = minuslog.lik,
start = list(p = hat))
summary(est)## Maximum likelihood estimation
##
## Call:
## stats4::mle(minuslogl = minuslog.lik, start = list(p = hat))
##
## Coefficients:
## Estimate Std. Error
## p 0.1357098 0.008564278
##
## -2 log L: 1269.921The maximum likelihood estimate is 0.136 and the standard error is 0.008564.