The following example is drawn from examples D and E in Section 4 in Chapter 4 of the Third Edition of Mathematical Statistics and Data Analysis by John A. Rice.

An Example

Say we want to estimate the total amount of money spent by all of the customers that enter a store in one day. We'll call this E(T). This amount is conditional on the amount of customers which enter the store on that day. This is N. Let's say that N is a random Poisson variable with a rate parameter \(\lambda\). Xi is the expected amount of money each of the N customers will spend. Each customer's spending is independent.


Finding the Expectation

Let's say that each customer is expected to spend the same amount E(X). Additionally, each customer's spending has the variance, Var(X). Because all Xis are independent, we know that

Using Adam's law and the fact above, we can conlude that


Finding the Variance

Because we know that

we also know that

With the knowledge that

we can deduce

Using Eve's Law we can finally conclude that

This app is designed to illustrate that these two conclusions are true with some simple simulations. The panels on the top of this page give different options for distributions of S, the amount each customer spends.


Theoretical Values

Comparison to Simulated Values


Theoretical Values

Comparison to Simulated Values


Theoretical Values

Comparison to Simulated Values

This App was created by Coco Kusiak to create a visualization of Adam's and Eve's laws designed to help students learning probability and its applications. Further information about these laws can be found at:

A special thanks to Professor Nicholas Horton and Professor Susan Wang.

Other resources include Mathematical Statisics and Data Analysis by John A. Rice and Introduction to Probability by Joseph K. Blitzstein and Jessica Hwang.

show with app
library(shiny)
library(mosaic)
library(lattice)
library(latticeExtra)


shinyServer(function(input, output) {
  
  
  FindingT <- function(Lamb, MIN, MAX){
    #Lamb <- input$Lamb
    #MIN <- input$MIN
    #MAX <- input$MAX
    N <- rpois(n = 1, lambda = Lamb)
    xSubi <- runif(n = N, min = MIN, max = MAX)
    sumXis <- sum(xSubi)
  }
  
  
  # display 10 rows initially
  tableData <- reactive({
    Simulation <- do(input$numsim)*FindingT(input$Lamb, input$MIN, input$MAX)
    ExpectationT <- mean(Simulation$FindingT)
    ExpectationN <- input$Lamb 
    ExpectationXi <- .5*(input$MIN + input$MAX)
    Theoretical_ExpectationT <- ExpectationN*ExpectationXi
    ForConfindence <- Simulation - Theoretical_ExpectationT
    confidence <- quantile(ForConfindence$FindingT, c(0.05, 0.95))
    confidence2 <- quantile(Simulation$FindingT, c(0.05, 0.95))
    c1 <- confidence2[1]
    c2 <- confidence2[2]
    Simulated_VarianceT <- var(Simulation$FindingT)
    VarianceXi <- (1/12)*(input$MAX - input$MIN)^2
    VarianceN <- input$Lamb
    Theoretical_VarianceT <- (ExpectationXi^2)*VarianceN + ExpectationN*VarianceXi
    list(Theoretical_ExpectationT=Theoretical_ExpectationT,
         ExpectationT=ExpectationT,
         Theoretical_VarianceT=Theoretical_VarianceT,
         Simulated_VarianceT=Simulated_VarianceT,
         ExpectationN=ExpectationN,
         VarianceN=VarianceN,
         ExpectationXi=ExpectationXi,
         VarianceXi=VarianceXi,
         Simulation=Simulation,
         Theoretical_ExpectationT=Theoretical_ExpectationT,
         c1=c1, 
         c2=c2)
  })
  output$d1 <- renderTable({
    data.frame(Theoretical = c("Expectation of N", "Variance of N", 
                               "Expectation of Xi", "Variance of Xi"), 
               Values = c(tableData()$ExpectationN, tableData()$VarianceN, 
                          tableData()$ExpectationXi, tableData()$VarianceXi))
  })
  output$d2 <- renderTable({
    data.frame(ValueofT = c("Theoretical", "Simulated"), 
               Expectation = c(tableData()$Theoretical_ExpectationT, tableData()$ExpectationT), 
               Variance = c(tableData()$Theoretical_VarianceT, tableData()$Simulated_VarianceT),
               StandardDeviation = c(sqrt(tableData()$Theoretical_VarianceT), sqrt(tableData()$Simulated_VarianceT)))
  })
  output$debug <- renderText({
    c("expectation:", tableData()$Theoretical_ExpectationT,
      "lower bound:", round(tableData()$c1, 3), "upper bound:",
      round(tableData()$c2, 3))
  })
  graphing <- function(points, conf, ex){
    cat("conf=", conf, "\n")
    a <- histogram(~points, xlab="E[T]",  main="Distribution of Simulated E[X]", col="slategray2")
    b <- latticeExtra :: layer(panel.abline(v=conf, col="red", lwd=3), data=list(conf=conf))
    c <- latticeExtra :: layer(panel.abline(v=ex, col="blue", lwd=3), data=list(ex=ex))
    d <- a + b + c
    return(d)
  }
  output$plotUnif <- renderPlot({
    datas <- tableData()$Simulation$FindingT
    lower <- tableData()$c1
    upper <- tableData()$c2
    conf <<- c(lower, upper)
    theoretical <- tableData()$Theoretical_ExpectationT
    graphing(points = datas, conf= c(lower, upper), ex= theoretical)
  })
  output$TEu1 <- renderUI({
    withMathJax(
      helpText("$$= \\left(\\frac{\\theta_{min}+\\theta_{max}}{2}\\right) \\lambda $$")
  )})
  output$TEu2 <- renderText({
    c("= ", tableData()$Theoretical_ExpectationT)
  })
  output$summation <- renderUI({
    withMathJax(
      helpText('$$T = \\sum_{i=1}^{N} X_i$$')
    )
  })
  output$expect1 <- renderUI({
    withMathJax(
      helpText('$$E[T|N] = NE[Xi]$$')
    )
  })
  output$expect2 <- renderUI({
    withMathJax(
      helpText("$$E(T) = E[NE(X)] = E(N)E(X)$$"),
      helpText("$$E(T) = E[NE[S]] = E[N]E(S)$$")
    )
  })
  output$var1 <- renderUI({
    withMathJax(
      helpText("$$Var(T|N) = N Var(X)$$")
    )
  })
  output$var2 <- renderUI({
    withMathJax(
      helpText("$$E[Var(T|N)] = E[N]Var(X)$$")
    )
  })
  output$var3 <- renderUI({
    withMathJax(
      helpText("$$E(T|N)=NE[X]$$")
    )
  })
  output$var4 <- renderUI({
    withMathJax(
      helpText("$$Var(E[T|N]) = E[X]^2 Var(N)$$")
    )
  })
  output$var5 <- renderUI({
    withMathJax(
      helpText("$$Var(T) = E[X]^2Var(N) + E[N]Var(X)$$")
    )
  })
  output$adams <- renderUI({
    withMathJax(
      helpText("$$E[T] = E[E[T|N]]$$")
    )
  })
  output$eves <- renderUI({
    withMathJax(
      helpText("$$Var(T) = E[Var(T|N)] + Var(E[T|N])$$")
    )
  })
  output$TVu1 <- renderUI({
    withMathJax(
      helpText(c("$$Var(T) = [E[X]]^2Var(N) + E[N]Var(X)$$")))
    })
  output$TVu2 <- renderUI({
    withMathJax(
      helpText(c("$$= \\left(\\frac{(\\theta_{min} + \\theta_{max})}{2}\\right)^2 \\lambda + 
                                      \\lambda  \\left( \\frac{(\\theta_{min} + \\theta_{max})^2}{12} \\right)^2$$")))
  })

  output$equations5LATEX <- renderUI({
    withMathJax(
      helpText(c("$$= \\left(\\frac{(", input$MIN,  "+", input$MAX, ")}{2}\\right)^2 \\lambda + 
                                      \\lambda  \\left( \\frac{(\\theta_{min} + \\theta_{max})^2}{12} \\right)^2$$")))
  })
  output$TVu3<- renderText({c("= (.5(", input$MIN, " + ", input$MAX, "))^2 x ", 
                                     input$lamb, " + ", input$lamb, " x (1/12(", input$MIN, " + ",
                                     input$MAX, ")^2)^2")})
  output$TVu4 <- renderText({c("= ", round(tableData()$Theoretical_VarianceT,2))})
  output$BelowPlot <- renderText("As shown above, the Theoretical E[T] (shown in blue), the total amount 
                                 spent at the store in one day, is within the bounds of a 95% confidence
                                 interval of the Simulated E[T] (shown in red).")
  output$BelowPlot2 <- renderText("As shown above, the Theoretical E[T] (shown in blue), the total amount 
                                 spent at the store in one day, is within the bounds of a 95% confidence
                                 interval of the Simulated E[T] (shown in red).")
  output$BelowPlot3 <- renderText("As shown above, the Theoretical E[T] (shown in blue), the total amount 
                                 spent at the store in one day, is within the bounds of a 95% confidence
                                 interval of the Simulated E[T] (shown in red).")
  FindingTNorm <- function(LambN, mean1, sd1){
    Lamb <- input$LambN
    mean1 <- input$mean1
    sd1 <- input$sd1
    N <- rpois(n = 1, lambda = LambN)
    xSubi <- rnorm(n = N,mean = mean1, sd = sd1)
    NormXsum <- sum(xSubi)
  }
  tableDataNorm <- reactive({
    Simu <- do(input$numsimN)*FindingTNorm(input$LambN, input$mean1, input$sd1)
    ExpectT <- mean(Simu$FindingTNorm)
    ExpectN <- input$LambN 
    ExpectXi <- input$mean1
    TExpectationT <- ExpectN*ExpectXi
    ForConfi <- Simu - TExpectationT
    confi<- quantile(ForConfi$FindingTNorm, c(0.05, 0.95))
    confi2 <- quantile(Simu$FindingTNorm, c(0.05, 0.95))
    c1N <- confi2[1]
    c2N <- confi2[2]
    SVarianceT <- var(Simu$FindingTNorm)
    VarXi <- (input$sd1)^2
    VarN <- input$LambN
    TVarianceT <- (ExpectXi^2)*VarN + ExpectN*VarXi
    list(TExpectationT=TExpectationT,
         ExpectT=ExpectT,
         TVarianceT=TVarianceT,
         SVarianceT=SVarianceT,
         ExpectN=ExpectN,
         VarN=VarN,
         ExpectXi=ExpectXi,
         VarXi=VarXi,
         Simu=Simu,
         TExpectationT=TExpectationT,
         c1N=c1N, 
         c2N=c2N)
  })
  
  output$d4 <- renderTable({
    data.frame(ValuesforT = c("Theoretical", "Simulated"), 
               Expectation = c(tableDataNorm()$TExpectationT, tableDataNorm()$ExpectT), 
               Variance = c(tableDataNorm()$TVarianceT, tableDataNorm()$SVarianceT), 
               StandardDeviation = c(sqrt(tableDataNorm()$TVarianceT), sqrt(tableDataNorm()$SVarianceT)))
  })
  output$d3 <- renderTable({
    data.frame(Theoretical = c("Expectation of N", "Variance of N", 
                               "Expectation of Xi", "Variance of Xi"), 
               Values = c(tableDataNorm()$ExpectN, tableDataNorm()$VarN, 
                          tableDataNorm()$ExpectXi, tableDataNorm()$VarXi))
  })
  output$plotNorm <- renderPlot({
    datas <- tableDataNorm()$Simu$FindingTNorm
    lower1 <- tableDataNorm()$c1N
    upper1 <- tableDataNorm()$c2N
    theory <- tableDataNorm()$TExpectationT
    graphing(points = datas, conf= c(lower1, upper1), ex= theory)
  })
  output$debug2 <- renderText({
    c("expectation:", tableDataNorm()$TExpectationT, "lower bound:", 
      round(tableDataNorm()$c1N, 3), "upper bound:", round(tableDataNorm$c2N, 3))
  })
  output$TEn1 <- renderUI({
    withMathJax(
      helpText("$$= \\mu\\lambda $$")
    )})
  output$TEn2 <- renderText({
    c("= ", tableDataNorm()$TExpectationT)
  })
  output$TVn1 <- renderUI({
    withMathJax(
      helpText(c("$$Var(T) = E[X]^2Var(N) + E[N]Var(X)$$")))
  })
  output$TVn2 <- renderUI({
    withMathJax(
      helpText(c("$$= \\mu^2\\lambda + \\lambda\\sigma^2$$")))
  })
  output$TVn3<- renderText({c("= (.5(", input$MIN, " + ", input$MAX, "))^2 x ", 
                              input$lamb, " + ", input$lamb, " x (1/12(", input$MIN, " + ",
                              input$MAX, ")^2)^2")})
  output$TVn4 <- renderText({c("= ", round(tableDataNorm()$TVarianceT,2))})
  FindingTGam <-
    function(LambG, shape1, rate1){
      LambG <- input$LambG
      shape1 <- input$shape1
      rate1 <- input$rate1
      NG <- rpois(n = 1, lambda = LambG)
      xSubiG <- rgamma(n = NG, rate = rate1, shape = shape1)
      GamXsum <- sum(xSubiG)
      return(GamXsum)
      
    }
  tableDataGamma <- reactive({
    SimuGam <- do(input$numsimG)*FindingTGam(input$LambG, input$shape1, input$rate1)
    ExT <- mean(SimuGam$FindingTGam)
    Mid <- median(SimuGam$FindingTGam)
    ExN <- input$LambG 
    ExXi <- input$shape1 / input$rate1
    TExT <- ExN*ExXi
    GamConf <- SimuGam - TExT
    Gconf2 <- quantile(SimuGam$FindingTGam, c(0.05, 0.95))
    c1G <- Gconf2[1]
    c2G <- Gconf2[2]
    SimVarG <- var(SimuGam$FindingTGam)
    VarGXi <- input$shape1/(input$rate1)^2
    VarGN <- input$LambG
    GamTVarianceT <- (ExXi^2)*VarGN + ExN*VarGXi
    list(TExT=TExT,
         ExT=ExT,
         GamTVarianceT=GamTVarianceT,
         SimVarG=SimVarG,
         ExN=ExN,
         VarGN=VarGN,
         ExXi=ExXi,
         VarGXi=VarGXi,
         SimuGam=SimuGam,
         TExT=TExT,
         c1G=c1G, 
         c2G=c2G)  
  })
  output$tableG1 <- renderTable({
    data.frame(Theoretical = c("Expectation of N", "Variance of N", 
                               "Expectation of Xi", "Variance of Xi"), 
               Values = c(tableDataGamma()$ExN, tableDataGamma()$VarGN, 
                          tableDataGamma()$ExXi, tableDataGamma()$VarGXi))
  })
  output$tableG2 <- renderTable({
    data.frame(ValueofT = c("Theoretical", "Simulated"), 
               Expectation = c(tableDataGamma()$TExT, tableDataGamma()$ExT), 
               Variance = c(tableDataGamma()$GamTVarianceT, tableDataGamma()$SimVarG),
               StandardDeviation = c(sqrt(tableDataGamma()$GamTVarianceT), sqrt(tableDataGamma()$SimVarG)))
  }
  )
  output$plotGam <- renderPlot({
    datas <- tableDataGamma()$SimuGam$FindingTGam
    lower1 <- tableDataGamma()$c1G
    upper1 <- tableDataGamma()$c2G
    theory <- tableDataGamma()$TExT
    graphing(points = datas, conf= c(lower1, upper1), ex= theory)
  })
  
  output$TEg1 <- renderUI({
    withMathJax(
      helpText("$$=\\left(\\frac{\\theta_{shape}}{\\theta_{rate}}\\right)  \\lambda $$")
    )})
  output$TEg2 <- renderText({
    c("= ", tableDataGamma()$TExT)
  })
  output$TVg1 <- renderUI({
    withMathJax(
      helpText(c("$$Var(T) = [E[X]]^2Var(N) + E[N]Var(X)$$")))
  })
  output$TVg2 <- renderUI({helpText(c("$$= \\left(\\frac{\\theta_{shape}}{\\theta_{rate}}\\right)^2 \\lambda + 
                                      \\lambda  \\left( \\frac{\\theta_{shape}}{\\theta_{rate}^2} \\right)$$"))
  })
  
  output$TVg3<- renderText({c("= (.5(", input$MIN, " + ", input$MAX, "))^2 x ", 
                              input$lamb, " + ", input$lamb, " x (1/12(", input$MIN, " + ",
                              input$MAX, ")^2)^2")})
  output$TVg4 <- renderText({c("= ", round(tableDataGamma()$GamTVarianceT,2))})
  
  
})
library(shiny)
library(ggplot2)



# Define UI for application that draws a histogram
shinyUI(navbarPage(
  title = "Adam's and Eve's laws",
  tabPanel("An Example",
           sidebarLayout(position ="right",
                         sidebarPanel(
                           withMathJax(),
                           h1("Adam's Law"),
                           p("Adam's Law or the Law of Total Expectation states that when 
                             given the coniditonal expectation of a random variable T which 
                             is conditioned on N, you can find the expected value of unconditional 
                             T with the following equation:"),
                           uiOutput('adams'),
                           hr(),
                           h1("Eve's Law"),
                           p("Eve's Law (EVVE's Law) or the Law of Total Variance is used to 
                             find the variance of T when it is conditional on N, 
                             it states that:"),
                           uiOutput('eves'),
                           hr(),
                           h4("Recall for this example, "),
                           p(strong("T"), " = total amount spent at the store"),
                           p(strong("N"), " = total number of customers that day"),
                           p(strong("Xi"), " = amount spent by the ith customer")
                           ),
                         mainPanel(
                           withMathJax(),
                           p("The following example is drawn from examples D and E in Section 4 in Chapter 4 of the Third Edition of", em(
                             "Mathematical Statistics and Data Analysis"), "by John A. Rice."),
                           h1("An Example"),
                           p("Say we want to estimate 
                             the total amount of money spent by all of the
                             customers that enter a store in one day. We'll call this E(T). This amount is conditional on 
                             the amount of customers which enter the store on that day. 
                             This is N. Let's say that N is a random Poisson variable 
                             with a rate parameter \\(\\lambda\\). Xi is the expected amount of money each of the 
                             N customers will spend. Each customer's spending is independent."
                           ),
                           uiOutput('summation'),
                           hr(),
                           h3("Finding the Expectation"),
                           p("Let's say that each customer is expected to spend the same amount 
                             E(X). Additionally, each customer's spending has the variance, Var(X). 
                             Because all Xis are independent, we know that "),

                           uiOutput('expect1'),
                           p("Using", strong("Adam's law"),"and the fact above, we can conlude that"),

                           uiOutput('expect2'),
                           hr(),
                           h3("Finding the Variance"),
                           
                           p("Because we know that"),
                           uiOutput('var1'),
                           
                           p("we also know that"),
                           uiOutput("var2"),

                           p("With the knowledge that"),
                           uiOutput('var3'),
                           
                           p("we can deduce"),
                           uiOutput('var4'),
                           
                           p("Using", strong("Eve's Law"),"we can finally conclude that"),
                           uiOutput('var5'),
                           
                           p("This app is designed to illustrate that these two conclusions are true with some 
                             simple simulations. The panels on the top of this page give different options for 
                             distributions of S, the amount each customer spends.")
                           )
                        )),
  
  tabPanel('Uniform',   
           sidebarPanel(
             withMathJax(),
             p("This example models Xi (the amount each customer spends) from a", strong("uniform distribution."), "
               Recall that N (the number of customers entering the store) comes from a",  
               strong("Poisson distribution."), " Enter parameter values below to define these distributions."),
             numericInput("MIN", "\\(\\theta_{min}\\): The Minimum Amount Spent at Store by Each Customer", "0"),
             numericInput("MAX", "\\(\\theta_{max}\\): The Maximum Amount Spent at Store by Each Customer", "100"),
             numericInput("Lamb", "\\(\\lambda\\): The Rate to Determine How Many Customers Enter the Store", 5),
             numericInput("numsim", "Number of simulations", 1000),
             hr(),
             h5("Theoretical Expectation of T"),
             p("$$E(T)$$"),
             p("$$= E[E(T|N)]$$"),
             p("$$= E(X) E(N)$$"), 
             uiOutput("TEu1"),
             textOutput("TEu2"),
             h5("Theoretical Variance of T"),
             uiOutput("TVu1"),
             uiOutput("TVu2"),
             textOutput("TVu4")
           ),
           mainPanel(  
                       plotOutput('plotUnif'),
                       textOutput('BelowPlot'),
                       hr(),
                       h4("Theoretical Values"),
                       tableOutput('d1'),
                       h4("Comparison to Simulated Values"),
                       tableOutput('d2'))),
  tabPanel('Normal',   
           sidebarPanel(
             p("This example models Xi (the amount each customer spends) from a", strong("normal distribution."), "
               Recall that N (the number of customers entering the store) comes from a",  
               strong("Poisson distribution."), " Enter parameter values below to define these distributions."),
             numericInput("mean1", "\\(\\mu\\): The Mean Amount Spent at Store by Each Customer", "50"),
             numericInput("sd1", "\\(\\sigma\\): The Standard Deviation Amount Spent at Store by Each Customer", "5"),
             numericInput("LambN", "\\(\\lambda\\): The Rate to Determine How Many Customers Enter the Store", 5),
             numericInput("numsimN", "Number of simulations", 1000),
             hr(),
             h5("Theoretical Expectation of T"),
             p("$$E(T)$$"),
             p("$$= E[E(T|N)]$$"),
             p("$$= E(X) E(N)$$"), 
             uiOutput("TEn1"),
             textOutput("TEn2"),
             h5("Theoretical Variance of T"),
             uiOutput("TVn1"),
             uiOutput("TVn2"),
             textOutput("TVn4")
           ),
           mainPanel(  
                       plotOutput('plotNorm'),
                       textOutput('BelowPlot2'),
                       hr(),
                       h4("Theoretical Values"),
                       tableOutput('d3'),
                       h4("Comparison to Simulated Values"),
                       tableOutput('d4'))),
  tabPanel('Gamma',   
           sidebarPanel(
             p("This example models Xi (the amount each customer spends) from a", strong("gamma distribution."), "
               Recall that N (the number of customers entering the store) comes from a",  
               strong("Poisson distribution."), " Enter parameter values below to define these distributions."),
             numericInput("shape1", "\\(\\theta_{shape}\\): Shape Parameter for theAmount Spent at Store by Each Customer", "500"),
             numericInput("rate1", "\\(\\theta_{rate}\\): Rate Parameter for the Amount Spent at Store by Each Customer", "10"),
             numericInput("LambG", "\\(\\lambda\\): The Rate to Determine How Many Customers Enter the Store", 5),
             numericInput("numsimG", "Number of simulations", 1000),
             hr(),
             h5("Theoretical Expectation of T"),
             p("$$E(T)$$"),
             p("$$= E[E(T|N)]$$"),
             p("$$= E(X) E(N)$$"), 
             uiOutput("TEg1"),
             textOutput("TEg2"),
             h5("Theoretical Variance of T"),
             uiOutput("TVg1"),
             uiOutput("TVg2"),
             textOutput("TVg4")
           ),
           mainPanel(  
                       plotOutput('plotGam'),
                       textOutput('BelowPlot3'),
                       hr(),
                       h4("Theoretical Values"),
                       tableOutput('tableG1'),
                       h4("Comparison to Simulated Values"),
                       tableOutput('tableG2')
                       
                      )),
  tabPanel('About',
           p("This App was created by Coco Kusiak to create a visualization 
             of Adam's and Eve's laws designed to help students learning probability and its applications. 
             Further information about these laws can be found at:"),
           wellPanel(
             helpText(   a("Adam's Law",     href="https://en.wikipedia.org/wiki/Law_of_total_expectation", target="_blank"), br(),
                         a("Eve's Law", href="https://en.wikipedia.org/wiki/Law_of_total_variance", target="_blank")
             )
           ),
           p("A special thanks to Professor Nicholas Horton and Professor Susan Wang."),
           p("Other resources include", em("Mathematical Statisics and Data Analysis"), 
             "by John A. Rice and ", em("Introduction to Probability"), "by Joseph K. 
                                        Blitzstein and Jessica Hwang."))
           
  
  
  
                           ))