library(shiny)
library(bslib)
library(dplyr)
library(tidyr)
library(ggplot2)

ui <- page_navbar(
  title = "Dice Uncertainty",
  theme = bs_theme(bootswatch = "flatly"),
  nav_panel(title = "Single",
            icon = icon("dice"),
            h5("Explore the uncertainty in estimating probability from various sample sizes, using the example of a 6-sided dice."),
            layout_sidebar(
              sidebar = sidebar(
                title = "Single Experiment",
                numericInput("num_rolls", "Number of Rolls: (1 - 10,000)", value = 1, min = 1, step = 1),
                actionButton("roll_button", "Roll")
              ),
              layout_columns(
                col_widths = c(4, 8),
                navset_card_tab(
                  title = "Dice Roll Outcomes",
                  nav_panel(title = "Table", tableOutput("roll_table")),
                  nav_panel(title = "Plot", plotOutput("bar_chart")),
                  full_screen = TRUE
                ),
                layout_columns(
                  col_widths = 12,
                  navset_card_tab(
                    title = "Summary of All Rolls",
                    nav_panel(
                      title = "Table",
                      tableOutput("summary_table"),
                      uiOutput("equation1")
                    ),
                    nav_panel(title = "Plot", plotOutput("bar_chart2")),
                    full_screen = TRUE
                  ),
                  card(
                    card_header("90% Confidence Interval, Normal Distribution"),
                    card_body(
                      uiOutput("equation2"),
                      uiOutput("equation3"),
                      uiOutput("equation4"),
                      uiOutput("equation5")
                    )
                  )
                )
              )
            )
  ),
  nav_panel(title = "Multiple",
            icon = icon("clipboard-list"),
            h5("Explore the uncertainty in estimating probability from various sample sizes, using the example of a 6-sided dice."),
            layout_sidebar(
              sidebar = sidebar(
                title = "Multiple Experiments",
                numericInput("num_rolls_right", "Number of rolls: (1 - 10,000)", value = 1, min = 1, step = 1),
                numericInput("num_experiments", "Number of experiments: (1 - 100)", value = 1, min = 1, step = 1),
                actionButton("run_button", "Run")
              ),
              layout_columns(
                col_widths = c(4, 8),
                card(
                  card_header("Estimated Probability of Rolling a 4 for Each Experiment"),
                  card_body(tableOutput("experiment_table"))
                ),
                layout_columns(
                  col_widths = 12,
                  card(
                    card_header("Probability of Rolling a 4 Across Experiments"),
                    card_body(plotOutput("line_plot")),
                    full_screen = TRUE
                  ),
                  card(
                    card_header("Summary of Rolling a 4 Across All Experiments"),
                    card_body(
                      tableOutput("experiment_summary_table"),
                      uiOutput("equation6"),
                      uiOutput("equation7"),
                      uiOutput("equation8")
                    )
                  )
                )
              )
            )
  ),
  nav_panel(
    title = "Multiple Comparison",
    icon = icon("layer-group"),
    h5("Explore the uncertainty in estimating probability from various sample sizes, using the example of a 6-sided dice."),
    layout_sidebar(
      sidebar = sidebar(
        title = "Multiple Experiments of Different Sample Sizes",
        numericInput("samplesize1", "Sample 1", min = 1, max = 10000, value = 12),
        numericInput("samplesize2", "Sample 2", min = 1, max = 10000, value = 100),
        numericInput("samplesize3", "Sample 3", min = 1, max = 10000, value = 1000),
        numericInput("samplesize4", "Sample 4", min = 1, max = 10000, value = 10000),
        h5(),
        numericInput("num_experiments2", "Number of experiments: (10 - 1000)", value = 20, min = 10, max = 1000, step = 1),
        actionButton("run_button2", "Run")
      ),
      layout_columns(
        col_widths = c(5, 7), 
        card(
          card_header("Estimated Probability of Rolling a 4 for Each Experiment"),
          navset_card_underline(
            title = NULL,
            nav_panel("Sample 1", tableOutput("experiment_table1")),
            nav_panel("Sample 2", tableOutput("experiment_table2")),
            nav_panel("Sample 3", tableOutput("experiment_table3")),
            nav_panel("Sample 4", tableOutput("experiment_table4")),
          ),
          full_screen = TRUE
        ),
        layout_columns(
          col_widths = 12,
          navset_card_tab(
            title = "Summary of Rolling a 4 Across All Experiments",
            nav_panel("Table", tableOutput("summary_table2")),
            nav_panel("Plot", plotOutput("bar_chart3")),
            full_screen = TRUE
          ),
          navset_card_tab(
            title = "Average Probability Distribution",
            nav_panel("Table", tableOutput("hist_table")),
            nav_panel("Plot", plotOutput("hist_plot")),
            nav_panel("Statistics", plotOutput("hist_stats_plot")),
            full_screen = TRUE
          )
        )
      )
    )
  )
)

server <- function(input, output, session) {
  
  # --- TAB 1 ---
  single_roll_results <- eventReactive(input$roll_button, {
    req(input$num_rolls)
    outcomes <- sample(1:6, input$num_rolls, replace = TRUE)
    data.frame(Roll = seq_len(input$num_rolls), Outcome = outcomes)
  }, ignoreNULL = FALSE)
  
  output$roll_table <- renderTable({ single_roll_results() })
  
  output$summary_table <- renderTable({
    df <- single_roll_results()
    req(df)
    summary <- data.frame(Outcome = 1:6, Count = as.vector(table(factor(df$Outcome, levels = 1:6))))
    summary$`Relative Frequency` <- formatC(summary$Count / sum(summary$Count), format="f", digits=3)
    summary
  })
  
  num_rolls_tab1 <- reactiveVal(NULL)
  observeEvent(input$roll_button, {
    num_rolls_tab1(input$num_rolls)
  })
  
  output$bar_chart <- renderPlot({
    df <- single_roll_results()
    req(df)
    counts <- table(factor(df$Outcome, levels = 1:6))
    par(mar = c(4.5, 4.5, 1.5, 1.5))
    max_ylim <- max(input$num_rolls*0.6, max(counts)*1.3)
    bar <- barplot(counts, names.arg = 1:6, xlab = "Outcome", ylab = "Count",
                   ylim = c(0, max_ylim), col = ifelse(1:6 == 4, "dodgerblue3", "salmon2"))
    abline(h = (num_rolls_tab1() / 6), col = "red", lty = 2)
    text(x = bar, y = counts, labels = counts, pos = 3)
  })
  
  output$bar_chart2 <- renderPlot({
    df <- single_roll_results()
    req(df)
    # n <- input$num_rolls
    n <- num_rolls_tab1()
    freqs <- table(factor(df$Outcome, levels = 1:6)) / n
    mu <- 1/6
    sd_val <- sqrt(mu * (1 - mu) / n)
    up <- mu + qnorm(0.95) * sd_val
    lo <- mu + qnorm(0.05) * sd_val
    
    par(mar = c(4.5, 4.5, 1.5, 1.5))
    max_ylim <- max(max(freqs, up)*1.5, 0.5)
    bar <- barplot(freqs, names.arg = 1:6, xlab = "Outcome", ylab = "Relative Frequency",
                   ylim = c(0, max_ylim), col = ifelse(1:6 == 4, "dodgerblue3", "salmon2"))
    abline(h = c(up, lo), col = "black", lty = 2, lwd = 2)
    
    # Restoring horizontal line text labels
    text(x = 0.5, y = up, labels = "90% Confidence", pos = 3, col = "black", cex = 0.8)
    text(x = bar, y = freqs, labels = formatC(freqs, format="f", digits=3), pos = 3)
  })
  
  # --- TAB 2 ---
  multi_experiment_results <- eventReactive(input$run_button, {
    req(input$num_rolls_right, input$num_experiments)
    n_rolls <- input$num_rolls_right
    n_exps <- input$num_experiments
    probs <- replicate(n_exps, sum(sample(1:6, n_rolls, replace = TRUE) == 4) / n_rolls)
    data.frame(
      Experiment = 1:n_exps,
      `Number of 4s Rolled` = probs * n_rolls,
      Probability = probs,
      `Square Error` = (probs - 1/6)^2,
      check.names = FALSE
    )
  })
  
  output$experiment_table <- renderTable({
    df <- multi_experiment_results()
    df %>% mutate(Probability = formatC(Probability, format="f", digits=3),
                  `Square Error` = formatC(`Square Error`, format="e", digits=2))
  })
  
  output$experiment_summary_table <- renderTable({
    df <- multi_experiment_results()
    data.frame(
      `Average Probability` = formatC(mean(df$Probability), format="f", digits=3),
      `True Probability` = formatC(1/6, format="f", digits=3),
      `Standard Error` = formatC(sd(df$Probability) / sqrt(nrow(df)), format="f", digits=4),
      check.names = FALSE
    )
  })
  
  output$line_plot <- renderPlot({
    df <- multi_experiment_results()
    
    n <- input$num_rolls_right
    mu <- 1/6
    sd_val <- sqrt(mu * (1 - mu) / n)
    up <- mu + qnorm(0.95) * sd_val
    lo <- mu + qnorm(0.05) * sd_val
    
    par(mar = c(4.5, 4.5, 1.5, 1.5))
    plot(
      df$Experiment, 
      df$Probability, 
      type = "o", 
      pch = 16, 
      cex = 1.75,
      xlab = "Experiment", 
      ylab = "Probability of Rolling a 4", 
      ylim = c(0, 0.5), 
      col = "dodgerblue3"
    )
    abline(h = 1/6, col = "red", lty = 2)
    abline(h = c(up, lo), col = "black", lty = 2, lwd = 2)
    
    # Restoring horizontal line text labels
    text(x = 2, y = up, labels = "90% Confidence", pos = 3, col = "black", cex = 0.8)
  })
  
  # --- TAB 3 ---
  comp_data <- eventReactive(input$run_button2, {
    sizes <- c(input$samplesize1, input$samplesize2, input$samplesize3, input$samplesize4)
    n_exps <- input$num_experiments2
    lapply(sizes, function(s) {
      probs <- replicate(n_exps, sum(sample(1:6, s, replace = TRUE) == 4) / s)
      data.frame(Experiment = 1:n_exps, Size = s, Probability = probs, 
                 `Square Error` = (probs - 1/6)^2, check.names = FALSE)
    })
  })
  
  format_comp_display <- function(df) {
    df %>% mutate(Probability = formatC(Probability, format="f", digits=3),
                  `Square Error` = formatC(`Square Error`, format="e", digits=2))
  }
  
  output$experiment_table1 <- renderTable({ format_comp_display(comp_data()[[1]]) })
  output$experiment_table2 <- renderTable({ format_comp_display(comp_data()[[2]]) })
  output$experiment_table3 <- renderTable({ format_comp_display(comp_data()[[3]]) })
  output$experiment_table4 <- renderTable({ format_comp_display(comp_data()[[4]]) })
  
  output$summary_table2 <- renderTable({
    bind_rows(comp_data()) %>%
      group_by(Size) %>%
      summarise(
        `Average Probability` = formatC(mean(Probability), format="f", digits=3),
        `True Probability` = formatC(1/6, format="f", digits=3),
        `Standard Error` = formatC(sd(Probability) / sqrt(n()), format="f", digits=4)
      ) %>% rename(`Sample Size` = Size)
  })
  
  output$hist_stats_plot <- renderPlot({
    df <- bind_rows(comp_data())
    
    skew_fun <- function(x) {
      x <- as.numeric(x)
      x <- x[is.finite(x)]
      n <- length(x)
      if (n < 3) return(NA_real_)
      s <- sd(x)
      if (!is.finite(s) || s == 0) return(NA_real_)
      sum(((x - mean(x)) / s)^3) * (n / ((n - 1) * (n - 2)))
    }
    
    stats_df <- df %>%
      group_by(Size) %>%
      summarise(
        Mean = mean(Probability, na.rm = TRUE),
        SD   = sd(Probability, na.rm = TRUE),
        Skew = abs(skew_fun(Probability)),   # abs(skew)
        .groups = "drop"
      ) %>%
      pivot_longer(
        cols = c("Mean", "SD", "Skew"),      # enforce column order before pivot
        names_to = "Metric",
        values_to = "Value"
      ) %>%
      mutate(
        Size = factor(Size, levels = c(input$samplesize1, input$samplesize2, input$samplesize3, input$samplesize4)),
        Metric = factor(Metric, levels = c("Mean", "SD", "Skew"))  # enforce bar/legend order
      )
    
    ggplot(stats_df, aes(x = Size, y = Value, fill = Metric)) +
      geom_col(position = position_dodge(width = 0.9), width = 0.8) +
      theme_minimal() +
      labs(
        x = "Sample Size",
        y = "Value",
        fill = "Statistic"
      ) +
      scale_fill_manual(values = c("Mean" = "dodgerblue3", "SD" = "firebrick3", "Skew" = "darkgreen"))
  })
  
  output$bar_chart3 <- renderPlot({
    df <- bind_rows(comp_data()) %>%
      group_by(Size) %>%
      summarise(Avg = mean(Probability), SE = sd(Probability) / sqrt(n())) %>%
      pivot_longer(cols = c(Avg, SE), names_to = "Metric", values_to = "Value")
    
    ggplot(df, aes(x = as.factor(Size), y = Value, fill = Metric)) +
      geom_bar(stat = "identity", position = "dodge") +
      scale_fill_manual(values = c("Avg" = "blue", "SE" = "red"), 
                        labels = c("Average Probability", "Standard Error")) +
      geom_hline(yintercept = 1/6, linetype = "dashed", color = "black") +
      theme_minimal() + labs(x = "Sample Size", y = "Value")
  })
  
  output$hist_table <- renderTable({
    wide <- bind_rows(comp_data()) %>%
      mutate(Size = paste("Size", Size)) %>%
      select(Experiment, Size, Probability) %>%
      pivot_wider(names_from = Size, values_from = Probability) %>%
      mutate(Experiment = as.character(Experiment))  # <-- key fix
    
    prob_cols <- setdiff(names(wide), "Experiment")
    
    # sample skewness; NA if not enough data or sd==0
    skew_fun <- function(x) {
      x <- as.numeric(x)
      x <- x[is.finite(x)]
      n <- length(x)
      if (n < 3) return(NA_real_)
      s <- sd(x)
      if (!is.finite(s) || s == 0) return(NA_real_)
      sum(((x - mean(x)) / s)^3) * (n / ((n - 1) * (n - 2)))
    }
    
    mean_row <- wide[1, , drop = FALSE]
    mean_row$Experiment <- "Mean"
    mean_row[prob_cols] <- lapply(prob_cols, \(nm) mean(wide[[nm]], na.rm = TRUE))
    
    sd_row <- wide[1, , drop = FALSE]
    sd_row$Experiment <- "Std Dev"
    sd_row[prob_cols] <- lapply(prob_cols, \(nm) sd(wide[[nm]], na.rm = TRUE))
    
    skew_row <- wide[1, , drop = FALSE]
    skew_row$Experiment <- "Skew"
    skew_row[prob_cols] <- lapply(prob_cols, \(nm) skew_fun(wide[[nm]]))
    
    out <- bind_rows(wide, mean_row, sd_row, skew_row)
    
    # optional formatting for display
    out %>%
      mutate(across(all_of(prob_cols), ~ formatC(as.numeric(.x), format = "f", digits = 4)))
  })
  
  output$hist_plot <- renderPlot({
    df <- bind_rows(comp_data())
    ggplot(df, aes(x = Probability, fill = as.factor(Size))) +
      geom_histogram(position = "dodge", bins = 25, alpha = 0.7) +
      scale_fill_manual(values = c("blue", "red", "green", "purple")) +
      theme_minimal() + labs(fill = "Sample Size", x = "Probability", y = "Frequency")
  })
  
  # --- EQUATIONS (UNCHANGED) ---
  output$equation1 <- renderUI({ withMathJax("$$\\text{Relative Frequency} = \\frac{\\text{Number of Occurrences}}{\\text{Number of Independent Trials}} = \\text{Estimated Probability}$$") })
  output$equation2 <- renderUI({ withMathJax("$$\\mu = \\frac{1}{6}$$") })
  output$equation3 <- renderUI({ withMathJax("$$\\sigma = \\sqrt{\\frac{\\mu(1 - \\mu)}{n_{rolls}}}$$") })
  output$equation4 <- renderUI({ withMathJax("$$\\text{Upper Limit} = \\mu + \\sigma \\cdot z_{0.95}$$") })
  output$equation5 <- renderUI({ withMathJax("$$\\text{Lower Limit} = \\mu + \\sigma \\cdot z_{0.05}$$") })
  output$equation6 <- renderUI({ withMathJax("$$\\text{Square Error} = (\\text{Observed Value} - \\text{Predicted Value})^2$$") })
  output$equation7 <- renderUI({ withMathJax("$$\\text{Standard Deviation (sample)} = s = \\sqrt{\\frac{\\sum_{i=1}^n (\\text{Observed}_i - \\text{Mean})^2}{n-1}}$$") })
  output$equation8 <- renderUI({ withMathJax("$$\\text{Standard Error} = SE = \\frac{s}{\\sqrt{n}}$$" ) })
}

shinyApp(ui, server)