例題: バイモーダルテスト

# Load the necessary libraries
#Original code by Valentina Tillman
library(mixtools)
library(mtbr)
library(diptest)

# Retrieve sample data
input_column <- commandArgs(trailingOnly = TRUE)
data <- mtb_get_column(input_column)
dip_test_result <- dip.test(data)

if (dip_test_result$p.value < 0.05) {
  # Fit a bimodal mixture model
  bimodal_fit <- normalmixEM(data, k = 2)
  
  # Manually extract parameter estimates and format them as a data frame
  bimodal_table <- data.frame(
    Mean = bimodal_fit$mu,
    Standard_Deviation = bimodal_fit$sigma,
    Proportion = bimodal_fit$lambda #tells you what % of the data is clustered around which mean. Also called lambda
  )
  
  # Define title and headers
  mytitle <- "Modeling a Bimodal Distribution"
  myheaders <- names(bimodal_table)
  
  # Add the table to the mtbr output
  mtb_add_table(columns = bimodal_table, headers = myheaders, title = mytitle)
  
  png("r_bimodal_image.png")
  plot(bimodal_fit, density = TRUE, which = 2)
  graphics.off()
  mtb_add_image("r_bimodal_image.png")
  
  # Now generate tolerance intervals using the parameters found using mixtools
  # Set the desired coverage level (e.g., 95%)
  coverage_level <- 0.95
  alpha <- 1 - coverage_level
  
  # Calculate the tolerance intervals for each component
  tolerance_intervals <- lapply(1:2, function(i) {
    mu <- bimodal_fit$mu[i]
    sigma <- bimodal_fit$sigma[i]
    n <- bimodal_fit$lambda[i]  # proportion of the component
    
    # Calculate the critical value for the normal distribution
    z <- qnorm(1 - alpha / (2 * n))
    
    # Calculate lower and upper bounds of the tolerance interval
    lower_bound <- mu - z * sigma
    upper_bound <- mu + z * sigma
    
    c(lower_bound, upper_bound)
  })
  
  # Show the tolerance intervals
  tolerance_intervals_df <- data.frame(
    Component = c("First Mode", "Second Mode"),
    Lower_Bound = sapply(tolerance_intervals, "[", 1),
    Upper_Bound = sapply(tolerance_intervals, "[", 2)
  )
  
  myheaders <- c("Component", "Lower Bound", "Upper Bound")
  mytitle <- "Tolerance Intervals for Bimodal Distribution"
  mtb_add_table(columns = tolerance_intervals_df, headers = myheaders, title = mytitle)
  
  
} else {
  mtb_add_message("This data is unimodal.")
}