Part 1

The faithful data set contains information about the Old Faithful geyser in Yellowstone National Park. Specifically, it records the length of each eruption (in minutes) and the waiting time until the next eruption (also in minutes).

In this practical, we’re interested in the waiting times between eruptions.

Import Necessary Libraries &Data

# Load the necessary libraries
library(datasets)
library(ggplot2)
library(mixtools)
## mixtools package, version 2.0.0, Released 2022-12-04
## This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772 and the Chan Zuckerberg Initiative: Essential Open Source Software for Science (Grant No. 2020-255193).
# Load the faithful dataset
data(faithful)

# Display the first few rows of the dataset
head(faithful)

This data is a 2-column data frame: β€’ eruptions: Length of eruptions (in mins) β€’ waiting: Time in between eruptions (in mins)

Question 1. Visualize the density of the variable of interest waiting using R.undefined

# Visualize the density of waiting times
ggplot(faithful, aes(x = waiting)) +
  geom_density(fill = "blue", alpha = 0.5) +
  labs(title = "Density Plot of Waiting Times",
       x = "Waiting Time (minutes)",
       y = "Density")

Question 2. When we plot the density of these waiting times, how many peaks do you notice? What does this suggest with respect to the waiting times?

Answer: We can see that the plot above is showing bimodal distribution and these peaks represents the 2 different modes of the activity of the geyser,The part in between 50 - 60 have a small peak which it indicates a short waiting time then there’s 80 have a big peak which indicates that it has a long waiting time

To investigate this further, we can decide to use a statistical technique called a Gaussian Mixture Model (GMM). This way of modeling data that comes from multiple groups (assuming 2 you see multiple peeks in the resulting plot in Question 1). The GMM works by assuming that each group of data follows a Gaussian (or normal) distribution. This is a common type of distribution that is determined by two parameters: the mean (the center of the distribution) and the variance (how spread out the distribution is).

Question 3. Estimate the parameters of your models. Hint: To estimate these parameters, we use the Expectation-Maximization (EM) algorithm. You can use the normalmixEM function from the mixtools package to estimate the parameters of the model. From the lecture, this is an iterative method that starts with sime initial guesses for the parameters, and then alternates between two steps:

  1. Expectation step (E-step): Based on the current parameter estimates, calculate the expected value of the log-likelihood. This is a measure of how likely the observed data are, given the parameters. In simpler terms, we use our current guess of the parameters to determine the β€œresponsibility” that each Gaussian distribution takes in explaining each observation in the data.

  2. Maximization step (M-step): Update the parameters to maximize the expected log-likelihood calculated in the E-step. In simpler terms, we update our guess of the parameters based on the results from the E-step.

We repeat these two steps until our estimates of the parameters stop changing significantly.

# Fit a Gaussian Mixture Model
set.seed(38)
gmm <- normalmixEM(faithful$waiting, k = 2)
## number of iterations= 38
# Display the estimated parameters
summary(gmm)
## summary of normalmixEM object:
##           comp 1    comp 2
## lambda  0.360887  0.639113
## mu     54.614901 80.091098
## sigma   5.871251  5.867711
## loglik at estimate:  -1034.002

Question 4. What does the output tell us? Explain each part (i.e., iterations,lambda, mu, etc.) of the output and its meaning.

  1. Iterations: The EM algorithm ran for 38 iterations before converging.
     
  2. Lambda (πœ†): The model found two components with weights 0.36 and
  0.64.This means the first component (Gaussian distribution) accounts 
  for 36% ofthe data, and the second accounts for 64%.
  
  3. Mu (πœ‡): The means of the two components are 54.6 and 80.1, suggesting
  the data clusters around these two values.
  
  4. Sigma (𝜎): The standard deviations are 5.9 and 5.9, indicating the
  spread around the means.
  
  5. Log-Likelihood: -1034.002 is the log-likelihood value of the final
  estimated model parameters.

Question 5. Visualize the estimation of your model on the same plot as the histogram of the data.

# Histogram and density of waiting times with GMM components
hist(faithful$waiting, breaks = 30, probability = TRUE, main = "Histogram and GMM of Waiting Times", xlab = "Waiting Time (minutes)")
lines(density(faithful$waiting), col = "blue", lwd = 2)
curve(gmm$lambda[1] * dnorm(x, mean = gmm$mu[1], sd = gmm$sigma[1]), col = "red", lwd = 2, add = TRUE)
curve(gmm$lambda[2] * dnorm(x, mean = gmm$mu[2], sd = gmm$sigma[2]), col = "green", lwd = 2, add = TRUE)
legend("topright", legend = c("Density", "Component 1", "Component 2"), col = c("blue", "red", "green"), lwd = 2)

Part 2

For mcdata.csv perform the hypothesis test that the mean is not equal to 454, i.e 𝐻0 : πœ‡ = 454 vs. 𝐻1 : πœ‡ β‰  454. Implement a Monte Carlo simulation to estimate the critical values.

Import Necessary Libraries & Data

# Load the necessary library
library(readr)

# Load the mcdata.csv dataset
mcdata <- read_csv("C:\\Users\\demon\\Downloads\\mcdata.csv")
## New names:
## Rows: 25 Columns: 2
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," dbl
## (2): ...1, mcdata
## β„Ή 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.
## β€’ `` -> `...1`
# Display the first few rows of the dataset
head(mcdata)

Question 1. Perform a t-test to test the null hypothesis.

# Perform a t-test
t_test <- t.test(mcdata$mcdata, mu = 454)
t_test
## 
##  One Sample t-test
## 
## data:  mcdata$mcdata
## t = -2.3584, df = 24, p-value = 0.02684
## alternative hypothesis: true mean is not equal to 454
## 95 percent confidence interval:
##  446.4995 453.5005
## sample estimates:
## mean of x 
##       450

Question 2. Use a Monte Carlo simulation to estimate the critical values. Generate a large number of random samples from a normal distribution with the same mean and standard deviation as our data, and calculate the t-statistic for each sample.

# Parameters for Monte Carlo simulation
set.seed(123)  # For reproducibility
num_simulations <- 10000  # Number of simulations
sample_size <- length(mcdata$mcdata)
sample_mean <- mean(mcdata$mcdata)
sample_sd <- sd(mcdata$mcdata)

# Generate random samples and calculate t-statistics
t_statistics <- replicate(num_simulations, {
  simulated_data <- rnorm(sample_size, mean = sample_mean, sd = sample_sd)
  t.test(simulated_data, mu = 454)$statistic
})

# Calculate the critical values at the 2.5% and 97.5% percentiles
critical_values <- quantile(t_statistics, c(0.025, 0.975))
critical_values
##       2.5%      97.5% 
## -4.6276364 -0.3972937

Question 3. From Q1, what does the t-test result indicate? Does your sample mean of your data different from the hypothesized mean? Would you reject the null hypothesis? Check the 95% confidence interval for the mean, is this another indication for your rejection of null hypothesis?

# Display the t-test result
t_test
## 
##  One Sample t-test
## 
## data:  mcdata$mcdata
## t = -2.3584, df = 24, p-value = 0.02684
## alternative hypothesis: true mean is not equal to 454
## 95 percent confidence interval:
##  446.4995 453.5005
## sample estimates:
## mean of x 
##       450
# Check the 95% confidence interval for the mean
t_test$conf.int
## [1] 446.4995 453.5005
## attr(,"conf.level")
## [1] 0.95

Question 4. The Monte Carlo simulation result gives you the critical values at the 2.5% and 97.5% percentiles of the distribution of t-statistics that you would expect if the null hypothesis were true. Does the t-statistic from your data fall outside this range as a result of your simulation? Does the Monte Carlo simulation suggests rejecting the null hypothesis?

# Check if the t-statistic from the data falls outside the critical values
t_statistic <- t_test$statistic
t_statistic
##         t 
## -2.358388
# Compare t-statistic with critical values
reject_null <- t_statistic < critical_values[1] || t_statistic > critical_values[2]
reject_null
## [1] FALSE