Stat Project_(Moutushi)

#1(A)

# Set parameters
N <- 1000
p <- 0.5
iterations <- 5000

# Initialize vector to store results
odds_ratios <- numeric(iterations)

# Simulate data and calculate odds ratios
for (i in 1:iterations) {
  # Generate sample of Bernoulli random variables
  sample_data <- rbinom(N, 1, p)
  
  # Calculate sample average
  sample_average <- mean(sample_data)
  
  # Calculate odds ratio
  odds_ratio <- sample_average / (1 - sample_average)
  
  # Store result
  odds_ratios[i] <- odds_ratio
}

# Plot histogram of simulated odds ratios
hist(odds_ratios, freq = FALSE, breaks = 30, main = "Histogram of Simulated Odds Ratios",
     xlab = "Odds Ratio", ylab = "Density")

# Calculate mean and variance of odds ratios
mean_odds_ratio <- mean(odds_ratios)
var_odds_ratio <- var(odds_ratios)

# Calculate standard error using Delta method
se_delta_method <- sqrt((1/(1-p)^3)*(p/N))

# Plot normal distribution with mean and standard error
curve(dnorm(x, mean = mean_odds_ratio, sd = se_delta_method), 
      col = "blue", lwd = 2, add = TRUE)

#1(B)

# Set parameters
N_values <- c(10, 30, 50, 100, 500)
p <- 0.5
iterations <- 5000

# Function to calculate standard error using Delta method
calculate_se_delta_method <- function(p, N) {
  return(sqrt((1/(1-p)^3)*(p/N)))
}

# Function to simulate data and plot histogram for given N
simulate_and_plot <- function(N) {
  # Initialize vector to store results
  odds_ratios <- numeric(iterations)
  
  # Simulate data and calculate odds ratios
  for (i in 1:iterations) {
    # Generate sample of Bernoulli random variables
    sample_data <- rbinom(N, 1, p)
    
    # Calculate sample average
    sample_average <- mean(sample_data)
    
    # Calculate odds ratio
    odds_ratio <- sample_average / (1 - sample_average)
    
    # Store result
    odds_ratios[i] <- odds_ratio
  }
  
  # Plot histogram of simulated odds ratios
  hist(odds_ratios, freq = FALSE, breaks = 30, 
       main = paste("Histogram of Simulated Odds Ratios (N =", N, ")"),
       xlab = "Odds Ratio", ylab = "Density")
  
  # Calculate mean and variance of odds ratios
  mean_odds_ratio <- mean(odds_ratios)
  var_odds_ratio <- var(odds_ratios)
  
  # Calculate standard error using Delta method
  se_delta_method <- calculate_se_delta_method(p, N)
  
  # Plot normal distribution with mean and standard error
  curve(dnorm(x, mean = mean_odds_ratio, sd = se_delta_method), 
        col = "blue", lwd = 2, add = TRUE)
  
  # Add legend
  legend("topright", legend = "Normal Curve", col = "blue", lwd = 2)
}

# Loop over different values of N
for (N in N_values) {
  simulate_and_plot(N)
}

#2(A)

# Function to compute g(p)
compute_g <- function(p) {
  return(p / (1 - p))
}

# Function to compute delta method approximation for expected value
compute_expected_value <- function(sample_mean) {
  return(sample_mean / (1 - sample_mean))
}

# Function to compute delta method approximation for variance
compute_variance <- function(sample_mean, N, p) {
  g_prime <- 1 / ((1 - p)^2)
  var_x_bar <- p * (1 - p) / N
  return((g_prime^2 * var_x_bar) / (p / (1 - p)))
}

# Set the sample sizes and true parameters
sample_sizes <- c(10, 30, 100, 1000)
true_parameters <- c(0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)

# Store results
results <- data.frame()

# Loop through each simulation configuration
for (N in sample_sizes) {
  for (p in true_parameters) {
    # Repeat 100 times for each configuration
    for (i in 1:100) {
      # Step 1: Sample Generation
      sample_data <- rbinom(N, 1, p)
      
      # Step 2: Compute Sample Mean
      sample_mean <- mean(sample_data)
      
      # Step 3: Calculate g(X_bar)
      g_x_bar <- compute_g(sample_mean)
      
      # Step 4: Compute Delta Method Approximation for Expected Value
      estimated_expected_value <- compute_expected_value(sample_mean)
      
      # Step 5: Compute Delta Method Approximation for Variance
      estimated_variance <- compute_variance(sample_mean, N, p)
      
      # Store results
      result <- data.frame(Sample_Size = N,
                           True_Parameter = p,
                           Estimated_Expected_Value = estimated_expected_value,
                           Estimated_Variance = estimated_variance)
      results <- rbind(results, result)
    }
  }
}

# Save results to CSV
write.csv(results, file = "delta_method_results_monday.csv", row.names = FALSE)

#2(A) (Another solution for finding Samples)

# Function to generate samples and compute delta method estimates
generate_samples <- function(sample_sizes, true_parameters) {
  # Function to compute g(p)
  compute_g <- function(p) {
    return(p / (1 - p))
  }
  
  # Function to compute delta method approximation for expected value
  compute_expected_value <- function(sample_mean) {
    return(sample_mean / (1 - sample_mean))
  }
  
  # Function to compute delta method approximation for variance
  compute_variance <- function(sample_mean, N, p) {
    g_prime <- 1 / ((1 - p)^2)
    var_x_bar <- p * (1 - p) / N
    return((g_prime^2 * var_x_bar) / (p / (1 - p)))
  }
  
  # Store results
  results <- list()
  
  # Loop through each simulation configuration
  for (N in sample_sizes) {
    for (p in true_parameters) {
      sample_data <- list()
      estimates <- list()
      
      # Repeat 100 times for each configuration
      for (i in 1:100) {
        # Step 1: Sample Generation
        sample_data[[i]] <- rbinom(N, 1, p)
        
        # Step 2: Compute Sample Mean
        sample_mean <- mean(sample_data[[i]])
        
        # Step 3: Calculate g(X_bar)
        g_x_bar <- compute_g(sample_mean)
        
        # Step 4: Compute Delta Method Approximation for Expected Value
        estimated_expected_value <- compute_expected_value(sample_mean)
        
        # Step 5: Compute Delta Method Approximation for Variance
        estimated_variance <- compute_variance(sample_mean, N, p)
        
        # Store estimates
        estimates[[i]] <- list(Estimated_Expected_Value = estimated_expected_value,
                               Estimated_Variance = estimated_variance)
      }
      
      # Store results for this configuration
      results[[paste("N_", N, "_p_", p, sep = "")]] <- list(Sample_Data = sample_data,
                                                            Estimates = estimates)
    }
  }
  
  return(results)
}

# Usage:
# Define sample sizes and true parameters
sample_sizes <- c(10, 30, 100, 1000)
true_parameters <- c(0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)

# Generate samples and compute estimates
samples <- generate_samples(sample_sizes, true_parameters)

# Access samples for a particular configuration
# For example, to access sample data for N = 10 and p = 0.5
sample_data_1 <- samples[["N_10_p_0.5"]][["Sample_Data"]]
sample_data_2 <- samples[["N_30_p_0.5"]][["Sample_Data"]]
sample_data_3<- samples[["N_100_p_0.5"]][["Sample_Data"]]
sample_data_4<- samples[["N_1000_p_0.5"]][["Sample_Data"]]
#sample_data
sample_data_2

[[1]]
 [1] 1 0 1 0 0 1 1 0 0 0 0 1 0 0 1 1 0 0 1 1 0 0 1 0 1 1 1 1 1 1

[[2]]
 [1] 1 0 1 1 0 0 0 0 0 1 0 1 0 1 0 0 1 1 0 1 1 1 0 1 0 1 1 0 1 1

[[3]]
 [1] 1 0 1 0 0 0 1 0 1 1 1 0 1 0 1 1 1 1 1 0 0 0 1 1 0 1 0 0 1 0

[[4]]
 [1] 1 0 0 0 1 1 0 1 0 0 0 0 1 1 0 1 1 1 0 0 0 0 1 1 1 1 0 0 1 0

[[5]]
 [1] 0 1 1 1 0 1 0 0 0 0 0 0 1 0 1 0 1 0 0 1 1 1 0 1 0 0 1 1 1 0

[[6]]
 [1] 0 0 1 0 1 1 1 1 0 1 1 0 0 0 1 0 0 0 1 1 0 1 1 1 0 1 1 1 1 1

[[7]]
 [1] 0 0 0 0 1 1 1 0 0 1 1 1 0 1 0 1 0 0 1 1 1 1 1 0 0 0 1 0 1 1

[[8]]
 [1] 0 1 0 1 0 0 0 0 1 0 0 1 1 0 1 1 1 1 0 1 0 0 0 1 0 0 0 0 1 1

[[9]]
 [1] 0 0 1 1 0 0 1 0 0 1 0 1 0 1 0 1 0 1 1 0 1 0 1 1 0 1 1 1 1 0

[[10]]
 [1] 1 0 0 1 1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 0 0 1 0 1 0 1 0 0 0

[[11]]
 [1] 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 1 1 1 0 1 1 0 0 1 0 1 0 1

[[12]]
 [1] 1 1 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 1 1 0 1 1 0 1 0 0

[[13]]
 [1] 0 1 0 0 0 0 0 0 1 0 1 1 0 1 1 1 0 0 1 1 1 0 0 1 0 0 1 0 0 1

[[14]]
 [1] 1 0 0 1 1 1 0 1 1 1 1 0 1 0 0 0 1 1 1 0 1 1 1 0 1 0 0 1 0 1

[[15]]
 [1] 0 1 0 0 0 0 1 0 0 1 1 1 0 1 1 0 0 1 1 1 1 1 1 0 0 0 0 0 1 0

[[16]]
 [1] 0 0 0 0 0 0 1 0 1 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 0 1 0 0

[[17]]
 [1] 1 0 0 1 1 1 0 0 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1 0 1 1 0 0 0 0

[[18]]
 [1] 1 0 1 1 0 1 0 0 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 1 1 0 1 0 0 0

[[19]]
 [1] 0 1 1 1 0 0 1 1 1 0 0 0 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0

[[20]]
 [1] 1 0 1 1 1 0 1 0 1 0 0 1 0 1 0 0 1 1 0 1 0 0 0 1 1 1 0 0 0 0

[[21]]
 [1] 1 0 0 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 0 1 0 0 0

[[22]]
 [1] 0 0 0 1 0 0 0 1 1 0 0 0 0 1 1 1 0 1 0 1 1 1 1 1 1 0 1 0 0 1

[[23]]
 [1] 1 0 0 1 0 0 0 1 1 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 1 1 0 0

[[24]]
 [1] 1 0 0 1 1 1 0 0 1 1 0 0 1 1 0 1 0 1 0 1 0 0 1 1 0 1 1 1 0 0

[[25]]
 [1] 1 0 1 1 0 1 0 0 0 0 1 1 0 0 1 1 0 1 1 1 0 1 0 0 1 1 0 1 1 1

[[26]]
 [1] 0 1 0 1 0 1 0 0 1 0 0 0 0 0 1 0 1 1 0 0 0 1 0 0 1 0 1 1 0 1

[[27]]
 [1] 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 1 1 1 1 1 1 0 0 1 0 1 1 0

[[28]]
 [1] 1 1 1 0 0 0 1 0 0 0 1 1 0 1 0 0 1 0 1 1 1 0 0 1 0 0 1 0 0 0

[[29]]
 [1] 0 1 0 0 0 0 0 1 0 1 0 1 1 0 0 0 1 0 1 1 1 1 0 0 1 1 1 1 0 1

[[30]]
 [1] 1 0 1 0 0 0 0 1 1 0 0 1 0 0 0 1 1 1 1 1 1 0 0 1 0 1 1 0 1 0

[[31]]
 [1] 0 0 0 1 0 1 0 0 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 0 1 0 0 0 1 0

[[32]]
 [1] 1 1 0 0 0 0 1 1 1 1 0 1 0 0 0 1 0 0 0 1 0 1 0 1 0 0 1 1 1 0

[[33]]
 [1] 0 1 0 1 0 0 1 1 1 0 0 0 0 1 1 1 0 1 0 1 0 0 1 0 0 0 0 1 1 0

[[34]]
 [1] 1 0 0 1 0 0 1 1 0 0 1 0 1 0 0 0 0 0 0 1 0 1 1 1 0 0 1 1 1 1

[[35]]
 [1] 1 0 0 0 1 1 1 0 0 1 0 1 1 0 1 1 0 1 0 1 0 0 0 0 0 1 0 0 1 1

[[36]]
 [1] 1 1 1 0 1 1 1 1 1 1 0 1 1 0 1 0 1 1 0 0 0 1 0 1 0 1 1 0 0 1

[[37]]
 [1] 0 1 1 0 1 0 1 1 1 0 0 1 0 0 1 1 1 1 1 1 1 1 0 0 1 0 0 0 0 0

[[38]]
 [1] 1 0 1 1 0 1 1 1 1 1 1 1 0 1 1 1 0 1 0 0 0 0 1 0 0 0 0 1 1 0

[[39]]
 [1] 1 0 0 1 0 1 1 1 0 1 1 0 1 0 1 0 1 0 0 0 1 0 1 1 1 0 0 0 1 0

[[40]]
 [1] 1 0 0 1 0 1 0 0 0 1 0 1 0 0 1 0 0 1 1 0 1 1 0 1 0 1 0 0 0 0

[[41]]
 [1] 1 0 0 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 0 0 0 1 1 0 1 0 0 0

[[42]]
 [1] 0 1 1 1 1 1 0 0 1 1 1 1 0 1 1 0 1 0 0 1 0 1 0 1 0 1 1 0 1 0

[[43]]
 [1] 0 0 1 1 1 1 0 1 0 1 0 0 1 1 0 1 0 0 1 1 0 0 0 0 1 1 1 1 0 1

[[44]]
 [1] 1 0 0 1 0 1 0 1 0 0 1 1 0 1 1 0 0 0 0 1 0 0 1 0 1 1 1 1 0 0

[[45]]
 [1] 1 0 0 0 1 0 0 0 1 1 1 1 0 0 0 1 0 1 0 0 0 0 0 1 1 1 1 0 1 0

[[46]]
 [1] 0 1 0 1 0 0 1 0 0 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0

[[47]]
 [1] 1 0 0 0 1 1 1 1 1 1 1 1 1 0 1 0 0 1 0 1 0 1 1 1 0 0 0 1 0 0

[[48]]
 [1] 1 1 1 0 0 0 0 0 1 1 1 1 1 1 1 0 1 0 1 1 0 1 1 0 0 0 1 0 1 0

[[49]]
 [1] 0 0 0 0 1 0 1 0 1 0 1 0 0 0 1 0 0 1 1 1 0 1 0 0 0 0 0 0 1 0

[[50]]
 [1] 0 1 0 0 1 0 0 1 1 1 1 1 1 0 1 0 1 1 0 1 0 0 1 0 1 1 1 1 1 0

[[51]]
 [1] 0 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 0 1 1 1 1 0 1 0 1

[[52]]
 [1] 1 0 1 1 1 0 0 1 0 0 1 0 1 0 0 0 0 1 0 0 0 1 1 1 0 0 0 1 0 1

[[53]]
 [1] 0 0 1 0 0 0 1 0 0 0 1 1 1 0 0 0 1 0 0 0 0 0 1 1 0 1 1 0 0 1

[[54]]
 [1] 1 1 1 0 0 1 1 0 0 1 0 0 1 0 0 1 0 1 1 0 1 1 1 1 0 1 0 1 0 0

[[55]]
 [1] 1 0 1 0 1 0 0 0 0 0 0 1 1 1 1 0 1 1 0 1 0 0 0 0 0 1 1 0 1 0

[[56]]
 [1] 1 1 1 0 0 1 1 1 0 1 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 1 0

[[57]]
 [1] 0 0 0 0 0 0 0 0 1 1 1 1 0 1 1 1 1 0 0 0 1 1 1 1 0 0 1 1 1 1

[[58]]
 [1] 0 0 0 1 1 0 1 0 0 1 0 0 1 1 1 0 0 0 0 0 1 1 1 0 1 0 0 1 0 0

[[59]]
 [1] 1 1 0 0 0 1 1 0 0 1 0 0 0 1 0 1 1 1 1 0 0 0 1 0 0 1 1 0 0 0

[[60]]
 [1] 1 1 1 1 0 1 0 0 1 1 0 0 1 1 0 1 0 1 1 0 0 0 0 1 1 0 1 0 0 0

[[61]]
 [1] 0 1 1 1 0 1 0 1 0 1 0 1 1 1 0 0 0 1 1 1 0 1 0 1 1 0 0 0 1 0

[[62]]
 [1] 1 1 1 0 1 0 0 1 0 1 0 1 1 1 0 0 0 0 0 1 1 1 0 1 1 1 0 1 0 0

[[63]]
 [1] 0 1 1 0 1 1 0 0 1 0 1 1 0 0 0 0 1 1 0 0 1 1 1 1 0 1 1 1 0 1

[[64]]
 [1] 0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0 1 1 1 1 1 0 1 0 1 1 0 0 0 0

[[65]]
 [1] 1 0 0 0 1 0 0 1 1 1 0 1 1 1 0 0 1 1 1 1 1 0 1 1 1 1 1 0 0 1

[[66]]
 [1] 0 1 0 0 1 1 0 1 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 0 0 0 0 0 0 1

[[67]]
 [1] 0 0 0 1 0 0 0 1 1 1 1 1 0 1 0 1 1 1 1 0 1 0 0 1 1 1 0 0 1 1

[[68]]
 [1] 0 1 0 0 1 1 0 1 1 0 0 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 1 0 1 0

[[69]]
 [1] 1 0 1 1 1 1 0 1 0 0 1 0 1 1 0 0 0 0 0 1 0 1 0 1 1 1 0 1 1 0

[[70]]
 [1] 1 1 0 1 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 1 0 0

[[71]]
 [1] 0 0 0 0 1 0 0 1 1 1 0 0 1 0 0 1 0 0 1 1 1 0 1 1 1 1 1 1 1 0

[[72]]
 [1] 0 1 1 1 0 1 0 1 1 0 0 0 1 0 1 1 1 0 0 0 1 0 1 0 0 0 0 1 0 1

[[73]]
 [1] 0 1 0 0 0 1 0 0 1 1 1 0 1 0 0 1 1 0 1 1 0 0 1 0 0 0 1 1 1 1

[[74]]
 [1] 0 1 1 0 0 1 1 1 0 1 0 1 0 1 0 1 1 1 1 0 0 0 0 0 1 0 0 1 0 1

[[75]]
 [1] 1 1 0 0 1 0 0 1 0 1 0 1 1 1 1 1 0 0 0 0 1 0 1 1 1 0 0 0 0 0

[[76]]
 [1] 1 0 1 0 1 1 0 1 0 0 0 1 1 1 1 0 0 1 0 0 1 0 1 0 1 1 1 0 0 1

[[77]]
 [1] 1 1 0 0 1 0 1 0 1 1 1 1 0 1 0 0 0 1 1 0 0 0 0 1 0 1 1 1 1 1

[[78]]
 [1] 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 1 0 0 0 0 0 1 1 1 0

[[79]]
 [1] 1 0 0 0 0 1 1 1 1 1 0 1 0 1 0 0 1 0 1 0 1 1 0 0 1 1 1 0 1 1

[[80]]
 [1] 0 0 0 1 0 1 1 1 1 0 0 0 1 0 1 1 0 0 0 0 0 0 1 0 1 0 0 1 0 1

[[81]]
 [1] 1 1 1 1 0 0 0 0 0 1 1 1 1 0 1 0 1 0 0 1 0 1 0 0 1 0 0 1 1 0

[[82]]
 [1] 1 1 1 1 0 1 0 1 0 0 1 1 1 1 0 1 0 1 0 0 0 0 0 1 1 1 0 0 1 0

[[83]]
 [1] 1 1 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 0 1 0 1 0 1 0 0 1

[[84]]
 [1] 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 1 1 1 1 1 0 0

[[85]]
 [1] 1 0 1 1 0 0 1 0 0 0 0 1 0 1 0 0 0 1 1 1 0 1 1 1 1 0 1 1 0 1

[[86]]
 [1] 0 0 1 1 1 0 1 1 0 0 1 0 1 0 0 0 1 0 0 1 1 0 0 1 0 0 1 1 1 1

[[87]]
 [1] 0 1 0 0 1 1 1 1 1 0 1 0 0 1 1 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1

[[88]]
 [1] 0 0 0 1 1 0 0 0 1 0 1 0 0 1 0 0 0 1 1 0 1 1 0 1 1 1 1 0 0 0

[[89]]
 [1] 1 0 1 0 0 1 0 1 0 0 1 1 0 1 1 0 1 0 0 1 0 0 0 0 0 0 1 0 1 0

[[90]]
 [1] 0 0 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 0 1 0 0 1 0

[[91]]
 [1] 0 1 1 0 0 1 1 1 0 1 0 0 1 0 0 1 0 0 0 1 0 1 1 0 1 1 1 1 1 1

[[92]]
 [1] 1 0 0 0 1 0 1 1 1 0 1 0 1 0 1 0 0 1 0 0 1 1 0 1 0 1 0 0 0 1

[[93]]
 [1] 0 1 1 0 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1 1 0 0 1 0 1 1 1 0 0 0

[[94]]
 [1] 1 0 0 0 0 1 0 1 0 0 1 0 1 1 0 1 1 1 1 1 0 1 0 1 0 1 0 0 1 0

[[95]]
 [1] 1 0 1 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 0 1 0 1 0 1 1 1 1

[[96]]
 [1] 1 1 0 1 0 1 0 0 0 0 0 0 1 1 0 1 0 0 1 1 1 1 1 0 1 0 1 0 0 0

[[97]]
 [1] 1 1 0 1 1 0 1 1 0 0 0 0 0 1 0 0 1 1 1 1 0 0 1 0 0 0 0 1 0 1

[[98]]
 [1] 1 1 1 1 1 0 1 1 1 0 0 0 0 1 1 0 1 1 0 1 1 1 0 1 1 0 1 0 0 1

[[99]]
 [1] 0 0 1 0 0 1 0 1 0 1 1 1 0 1 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0

[[100]]
 [1] 0 1 0 0 0 1 1 1 1 1 1 0 0 0 1 0 1 1 0 1 0 1 1 1 1 0 1 0 1 0