Skip to content

Instantly share code, notes, and snippets.

@calpolystat

calpolystat/#MLE_Binomial.txt

Last active July 29, 2022 09:14
Show Gist options
  • Save calpolystat/8ace862e9e43f8e29d43 to your computer and use it in GitHub Desktop.
Save calpolystat/8ace862e9e43f8e29d43 to your computer and use it in GitHub Desktop.
Download ZIP
MLE for the Binomial distribution: Shiny app at http://www.statistics.calpoly.edu/shiny
Raw
#MLE_Binomial.txt
MLE Binomial Shiny App
Base R code created by Gail Potter
Shiny app files created by Gail Potter
Cal Poly Statistics Dept Shiny Series
http://statistics.calpoly.edu/shiny
Raw
DESCRIPTION
Title: MLE Binomial Distribution
Author: Gail Potter
AuthorUrl: http://www.gailpotter.org
License: MIT
DisplayMode: Normal
Tags: MLE Binomial Distribution
Type: Shiny
Raw
LICENSE
The MIT License (MIT)
Copyright (c) 2015 Gail Potter
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
Raw
server.R
# ------------------
# App Title: MLE Binomial
# Author: Gail Potter
# ------------------
shinyServer(
function(input, output, session) {
observe({
if (input$n > 0 & input$n %% 1 ==0& is.numeric(input$n) & !is.na(input$n))
return()
showshinyalert(session, "shinyalert1",
paste("Please enter a positive integer for n:"))
})
observe({
if (input$prob >= 0 & input$prob <=1 & is.numeric(input$prob))
return()
showshinyalert(session, "shinyalert2",
paste("Please enter a number between 0 and 1:"))
})
observe({
if (input$n2 > 0 & input$n2 %% 1 ==0& is.numeric(input$n2) & !is.na(input$n2))
return()
showshinyalert(session, "shinyalert3",
paste("Please enter a positive integer for n:"))
})
observe({
if (input$x <= input$n2 & input$x %% 1 ==0& is.numeric(input$x) &
!is.na(input$n2) & input$x >=0)
return()
showshinyalert(session, "shinyalert4",
paste("Please enter an integer between 0 and ", input$n2, ":", sep=""))
})
observe({
if (input$n3 > 0 & input$n3 %% 1 ==0& is.numeric(input$n3) & !is.na(input$n3))
return()
showshinyalert(session, "shinyalert5",
paste("Please enter a positive integer for n:"))
})
observe({
if (input$x2 <= input$n3 & input$x2 %% 1 ==0& is.numeric(input$x2) &
!is.na(input$n3) & input$x2 >=0)
return()
showshinyalert(session, "shinyalert6",
paste("Please enter an integer between 0 and ", input$n3, ":", sep=""))
})
output$pmf <- renderPlot({
# draw.plotmath.cell(expression (bgroup("(", atop(x,y), ")")))
if (!(input$n <= 0 | input$n %% 1 !=0)) n=input$n else n=10
if (input$prob >= 0 & input$prob <=1 & is.numeric(input$prob)) prob = input$prob else prob=0.5
plot(0:n, dbinom(0:n, size=n, prob=prob), col="blue",pch=16,
xlab="X (number of successes)", ylab="", main=paste("Binomial Probability Mass Function, P(X=x), when p=", prob, "n=",n))
segments(x0=0:n, y0=rep(0, (n+1)), x1=0:n, y1=dbinom(0:n, size=n, prob=prob),lwd=2,col="blue")
abline(h=0)
})
output$lik = renderPlot({
p=seq(0,1,by=0.0001)
if (input$x <= input$n2 & input$x %% 1 ==0& is.numeric(input$x) &
!is.na(input$n2) & input$x >=0) x=input$x else x=0
if (input$n2 > 0 & input$n2 %% 1 ==0&
is.numeric(input$n2) & !is.na(input$n2)) n=input$n2 else n2=10
plot(p, dbinom(x=x, size=n, prob=p), type="l", ylab="", main=
paste("Binomial likelihood function, L(p|X=", x, ", ", "n=", n, ")", sep="") )
abline(v=x/n, col="red", lty=2)
})
output$loglik = renderPlot({
p=seq(0,1,by=0.0001)
if (input$x2 <= input$n3 & input$x2 %% 1 ==0& is.numeric(input$x2) &
!is.na(input$n3) & input$x2 >=0) x=input$x2 else x=0
if (input$n3 > 0 & input$n3 %% 1 ==0&
is.numeric(input$n3) & !is.na(input$n3)) n=input$n3 else n=10
par(mfrow=c(1,2))
plot(p, dbinom(x=x, size=n, prob=p), type="l", ylab="", main=
paste("Likelihood function, \n L(p|X=", x, ", ", "n=", n, ")", sep="") )
abline(v=x/n, col="red", lty=2)
plot(p, log(dbinom(x=x, size=n, prob=p)), type="l", ylab="", main=
paste("Log likelihood function, \n log(L(p|X=", x, ", ", "n=", n, "))", sep="") )
abline(v=x/n, col="red", lty=2)
})
}
)
Raw
ui.R
# ------------------
# App Title: MLE Binomial
# Author: Gail Potter
# ------------------
# A couple IMPORTANT notes on mathJax:
# 1) withMathJax right now is only working with the developer version of Shiny,
# so that is needed in order for the math expressions to display properly.
# 2) the math expressions display in a webpage but not in the RStudio that opens
# up when testing the app. click "Open in browser" for proper display.
if (!require("devtools"))
install.packages("devtools")
if (!require("shinysky")) devtools::install_github("ShinySky","AnalytixWare")
library(shinysky)
shinyUI(navbarPage("Maximum Likelihood Estimation for the Binomial distribution",
tabPanel("Probability Mass Function",
tags$head(tags$link(rel = "icon", type = "image/x-icon", href =
"https://webresource.its.calpoly.edu/cpwebtemplate/5.0.1/common/images_html/favicon.ico")),
fluidRow(
column(3,
wellPanel(
helpText("Parameters for the binomial probability mass function"),
shinyalert("shinyalert1", TRUE, auto.close.after=10),
numericInput("n", label="Number of trials", value=10),
shinyalert("shinyalert2", TRUE, auto.close.after=10),
numericInput("prob", label="Probability of success", value=0.5, min=0, max=1),
br(),br(),br(),
div("Shiny app by",
a(href="http://www.gailpotter.org",target="_blank",
"Gail Potter"),align="right", style = "font-size: 8pt"),
div("Base R code by",
a(href="http://www.gailpotter.org",target="_blank",
"Gail Potter"),align="right", style = "font-size: 8pt"),
div("Shiny source files:",
a(href="https://gist.github.com/calpolystat/8ace862e9e43f8e29d43",
target="_blank","GitHub Gist"),align="right", style = "font-size: 8pt"),
div(a(href="http://www.statistics.calpoly.edu/shiny",target="_blank",
"Cal Poly Statistics Dept Shiny Series"),align="right", style = "font-size: 8pt")
)),
column(9,
wellPanel(h3("Probability Mass Function"),
withMathJax(p("You have learned about the ", em("probability mass function"), "(PMF) for the binomial random variable.
This is a function which has two parameters, n (number of trials) and p (probability of success), which determine its shape.
The function takes as input the number of successes, and gives as output the probability of that many successes in n trials.
The equation for the probability mass function is $$P(X=x) ={n \\choose x} p^x (1-p)^{n-x}$$")),
plotOutput("pmf"),
br(),
p("In an experiment, you usually don't know which of these possible PMFs is the truth, and
you observe a single value of x, the number of successes.
The goal is to estimate p based on your observation, x.
This motivates the ", em("likelihood function."), " In the likelihood function, the functional form
is the same, but we treat p as variable and x, as fixed. The number of trials, n, is also fixed (by the experimental design).")
))
)),
tabPanel("Likelihood Function",
fluidRow(
column(3,
wellPanel(
helpText("Specifications for the likelihood function:"),
shinyalert("shinyalert3", TRUE, auto.close.after=10),
numericInput("n2", label="Number of trials", value=10),
shinyalert("shinyalert4", TRUE, auto.close.after=10),
numericInput("x", label="Number of successes", value=0),
br(), br(), br(),
div("Shiny app by",
a(href="http://www.gailpotter.org",target="_blank",
"Gail Potter"),align="right", style = "font-size: 8pt"),
div("Base R code by",
a(href="http://www.gailpotter.org",target="_blank",
"Gail Potter"),align="right", style = "font-size: 8pt"),
div("Shiny source files:",
a(href="https://gist.github.com/calpolystat/8ace862e9e43f8e29d43",
target="_blank","GitHub Gist"),align="right", style = "font-size: 8pt"),
div(a(href="http://www.statistics.calpoly.edu/shiny",target="_blank",
"Cal Poly Statistics Dept Shiny Series"),align="right", style = "font-size: 8pt")
)),
column(9,
wellPanel(
h3("Likelihood function"),
withMathJax(p("Since the likelihood function treats p as a variable and x as fixed, it is written L(p|x, n).
The likelihood function for this example is $$L(p|x, n) ={n \\choose x} p^x (1-p)^{n-x} $$")),
br(),
plotOutput("lik"),
br(),
p("For different values of x and n, determine the value of p where the likelihood function achieves its max.
This value of p is called the Maximum Likelihood Estimate (MLE) for p.
Can you derive the formula for this estimator of p, in terms of x and n?")
)))
),
tabPanel("Log Likelihood Function",
fluidRow(
column(3,
wellPanel(
helpText("Specifications for the log likelihood function:"),
shinyalert("shinyalert5", TRUE, auto.close.after=10),
numericInput("n3", label="Number of trials", value=10),
shinyalert("shinyalert6", TRUE, auto.close.after=10),
numericInput("x2", label="Number of successes", value=0),
br(), br(), br(),
div("Shiny app by",
a(href="http://www.gailpotter.org",target="_blank",
"Gail Potter"),align="right", style = "font-size: 8pt"),
div("Base R code by",
a(href="http://www.gailpotter.org",target="_blank",
"Gail Potter"),align="right", style = "font-size: 8pt"),
div("Shiny source files:",
a(href="https://gist.github.com/calpolystat/8ace862e9e43f8e29d43",
target="_blank","GitHub Gist"),align="right", style = "font-size: 8pt"),
div(a(href="http://www.statistics.calpoly.edu/shiny",target="_blank",
"Cal Poly Statistics Dept Shiny Series"),align="right", style = "font-size: 8pt")
)),
column(9,
wellPanel(
h3("Log Likelihood function"),
withMathJax(p("In practice, it is usually easier to find the MLE by maximizing the ", em("log likelihood function"),
" instead of the likelihood function. Since the log transformation is monotone,
the two functions achieve their maximum in the same place. The log likelihood function for this example is
$$ \\log (L(p|x, n)) = \\log \\Big( {n \\choose x} p^x (1-p)^{n-x} \\Big) $$"
)),
plotOutput("loglik"),
br(),
p("We have introduced the concept of maximum likelihood in the context of estimating a binomial proportion, but the concept
of maximum likelihood is very general. Maximum likelihood is used to estimate parameters for a wide variety of distributions.")
)),
tags$style(type="text/css",
".shiny-output-error { visibility: hidden; }",
".shiny-output-error:before { visibility: hidden; }"
))
)
))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment