---
title: "Statistical distributions"
runtime: shiny
output: html_document
layout: raw
---

```{r, echo=FALSE}
fluidRow(
column(2,
  ## panel selection
  radioButtons(
      "distr", "Distribution",
      c(Bernoulli = "Bernoulli",
        Binomial = "Binomial",
        Poisson = "Poisson",
        Normal = "Normal",
        Lognormal = "Lognormal",
        Uniform = "Uniform",
        Beta = "Beta",
        Gamma = "Gamma"))),
column(5,
  ## these are shown for all panels
  sliderInput("n", label = "Sample size",
              min = 10, max = 1000, value = 100, step = 10),
  sliderInput("seed", label = "Random seed",
              min = 0, max = 100, value = 0, step = 10)),
column(5,
  ## Bernoulli
  conditionalPanel(
    condition = "input.distr == 'Bernoulli'",
      sliderInput("p", label = "Probability",
              min = 0, max = 1, value = 0.3, step = 0.05)),
  ## Binomial
  conditionalPanel(
    condition = "input.distr == 'Binomial'",
      sliderInput("p", label = "Probability",
              min = 0, max = 1, value = 0.3, step = 0.05),
      sliderInput("size", label = "Size",
              min = 1, max = 1000, value = 10, step = 50)),
  ## Poisson
  conditionalPanel(
    condition = "input.distr == 'Poisson'",
      sliderInput("lambda", label = "Mean/Rate",
              min = 0, max = 100, value = 5, step = 5)),
  ## Normal
  conditionalPanel(
    condition = "input.distr == 'Normal'",
      sliderInput("mu", label = "Mean",
              min = -10, max = 10, value = 0, step = 1),
      sliderInput("var", label = "Variance",
              min = 0.001, max = 10, value = 1, step = 0.5)),
  ## Logormal
  conditionalPanel(
    condition = "input.distr == 'Lognormal'",
      sliderInput("mu", label = "Mean",
              min = -10, max = 10, value = -1, step = 1),
      sliderInput("var", label = "Variance",
              min = -10, max = 10, value = 1, step = 1)),
  ## Uniform
  conditionalPanel(
    condition = "input.distr == 'Uniform'",
      sliderInput("a", label = "Minimum",
              min = -10, max = 10, value = -1, step = 0.5),
      sliderInput("b", label = "Maximum",
              min = -10, max = 10, value = 1, step = 0.5)),
  ## Beta
  conditionalPanel(
    condition = "input.distr == 'Beta'",
      sliderInput("shape1", label = "Shape 2",
              min = 0, max = 10, value = 1, step = 0.5),
      sliderInput("shape2", label = "Shape 1",
              min = 0, max = 10, value = 1, step = 0.5)),
  ## Gamma
  conditionalPanel(
    condition = "input.distr == 'Gamma'",
      sliderInput("shape", label = "Shape",
              min = 0.001, max = 10, value = 1, step = 0.5),
      sliderInput("rate", label = "Rate",
              min = 0.001, max = 10, value = 1, step = 0.5))))
fluidRow(renderPlot({
    par(las = 1)
    set.seed(input$seed)
    if (input$distr == "Uniform" && input$b < input$a)
        stop("Maximum must be greater than Minimum")
    y <- switch(input$distr,
        "Bernoulli" = rbinom(1000, 1, input$p),
        "Binomial" = rbinom(1000, input$size, input$p),
        "Poisson" = rpois(1000, input$lambda),
        "Normal" = rnorm(1000, input$mu, sqrt(input$var)),
        "Lognormal" = rlnorm(1000, input$mu, sqrt(input$var)),
        "Uniform" = runif(1000, input$a, input$b),
        "Beta" = rbeta(1000, input$shape1, input$shape2),
        "Gamma" = rgamma(1000, input$shape, input$rate))
    #plot(pt, L, type = "l", col="#3498db",
    #    ylab = "Likelihood", xlab="p",
    #    sub=paste0("Mean = ", round(mean(y[1:input$n]), 2), " (", 
    #        sum(1-y[1:input$n]), " 0s & ", sum(y[1:input$n]), " 1s)"),
    #    main = paste("Estimate =", round(pt[which.max(L)], 2)))
    #abline(v = input$p, lwd = 2, col = "#c7254e")
    #abline(v = pt[which.max(L)], lwd = 2, col = "#18bc9c")
    yy <- y[1:input$n]
    x <- switch(input$distr,
        "Bernoulli" = c(0,1),
        "Binomial" = seq(0, max(yy)+1, by = 1),
        "Poisson" = seq(0, max(yy)+1, by = 1),
        "Normal" = seq(min(yy)-1, max(yy)+1, length.out = 1000),
        "Lognormal" = seq(0.0001, max(yy)+1, length.out = 1000),
        "Uniform" = seq(input$a+0.0001, input$b-0.0001, length.out = 1000),
        "Beta" = seq(0.0001, 0.9999, length.out = 1000),
        "Gamma" = seq(0.0001, max(yy), length.out = 1000))
    d <- switch(input$distr,
        "Bernoulli" = dbinom(x, 1, input$p),
        "Binomial" = dbinom(x, input$size, input$p),
        "Poisson" = dpois(x, input$lambda),
        "Normal" = dnorm(x, input$mu, sqrt(input$var)),
        "Lognormal" = dlnorm(x, input$mu, sqrt(input$var)),
        "Uniform" = dunif(x, input$a, input$b),
        "Beta" = dbeta(x, input$shape1, input$shape2),
        "Gamma" = dgamma(x, input$shape, input$rate))
    xlab <- "x"
    ylab <- "Density"
    main <- paste0(input$distr, " distribution (n = ", input$n, ")")
    if (input$distr %in% c("Bernoulli", "Binomial", "Poisson")) {
        tmp <- table(yy) / input$n
        plot(tmp, ylim=c(0, max(tmp, d)), 
            ylab = ylab, xlab = xlab, main = main,
            col = "#cccccc", lwd = 10)
        points(x, d, pch = 21, col = "#c7254e", type = "b", 
            lty = 2, cex = 2)        
    } else {
        tmp <- hist(yy, plot = FALSE)
        hist(yy, freq = FALSE, ylim=c(0, max(tmp$density, d)), 
            ylab = ylab, xlab = xlab, main = main,
            col = "#ecf0f1", border = "#cccccc")
        lines(x, d, lwd = 2, col = "#c7254e")
    }
}))
```