Random Variable Generation: Shiny app at http://www.statistics.calpoly.edu/shiny

#RandVarGen.txt

Random Variable Generation Shiny App

Base R code created by Peter Chi
Shiny app files created by Peter Chi

Cal Poly Statistics Dept Shiny Series 
http://statistics.calpoly.edu/shiny

DESCRIPTION

Title: Random Variable Generation
Author: Peter Chi
AuthorUrl: http://statweb.calpoly.edu/~pchi
License: MIT
DisplayMode: Normal
Tags: Random variable generation, accept-reject, probability integral transform
Type: Shiny

LICENSE

The MIT License (MIT)
 
Copyright (c) 2015 Peter Chi
 
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.

rvg.R

#############################################
# Probability Integral Transform
#

# function to initialize dataframe, for either example
init.all<-function(){
  return(data.frame(NULL))
}

# Exponential
do.one.exp<-function(lambda,all.out){
  n.samples<-dim(all.out)[1]
  U<-runif(1,0,1)
  X<-(-log(1-U)/lambda)
  if(n.samples>0){
    height.X<-sum(round(all.out$X[1:n.samples],digits=2)==round(X,digits=2))
    height.U<-sum(round(all.out$U[1:n.samples],digits=2)==round(U,digits=2))
  } else {
    height.X<-0
    height.U<-0
  }
  return(data.frame(U=U,X=X,height.U=height.U,height.X=height.X))
}


plot.exp<-function(lambda,history){
  pdf.max<-lambda
  x.max<-max(c(1,qexp(0.95,lambda)))
  x<-seq(0,x.max,by=0.01)
  plot(dexp(x,lambda)~x,type='l',xlab="x",ylab="density",ylim=c(0,1.2*pdf.max))
  lines(c(0,0),c(0,pdf.max))

  n.points<-length(history$U)
  if(n.points>0){
#    for(i in 1:length(history$X)){
#       height[i]<-sum(round(history$X[1:i],digits=2)==round(history$X[i],digits=2))-1  # can i do this without a loop?
#    }
       color<-rep("black",n.points)
       color[n.points]<-"green"      # make the last one red
       cex<-rep(1.5,n.points) 
       cex[n.points]<-3
       hght.X<-history$height.X/(ifelse(max(history$height.X)<(10*pdf.max),10,(max(history$height.X)/pdf.max)))
       
       points(history$X,hght.X,pch="*",cex=cex,col=color)
   }
}


add.exp<-function(exp.all.out,lambda,num){
  for(i in 1:num){
    exp.all.out<-rbind(exp.all.out,do.one.exp(lambda,exp.all.out))    # need to loop this so that heights get added sequentially
  }
  return(exp.all.out)
#  return(rbind(exp.all.out,rdply(num,do.one.exp(lambda,exp.all.out))[,-1]))  # [,-1] is to get rid of the first column (which is indices)
}


# Linear Example
do.one.linear<-function(all.out){
  n.samples<-dim(all.out)[1]
  U<-runif(1,0,1)
  X<-(4*sqrt(U))
  if(n.samples>0){
    height.X<-sum(round(all.out$X[1:n.samples],digits=2)==round(X,digits=2))
    height.U<-sum(round(all.out$U[1:n.samples],digits=2)==round(U,digits=2))
  } else {
    height.X<-0
    height.U<-0
  }
  return(data.frame(U=U,X=X,height.U=height.U,height.X=height.X))
}



plot.linear<-function(history){
  pdf.max<-(1/2)
  x<-seq(0,4,by=0.01)
  plot((x/8)~x,type='l',xlab="x",ylab="density",ylim=c(0,1.2*pdf.max))
  lines(c(4,4),c(0,pdf.max))
  
  n.points<-length(history$U)
  if(n.points>0){
    #    for(i in 1:length(history$X)){
    #       height[i]<-sum(round(history$X[1:i],digits=2)==round(history$X[i],digits=2))-1  # can i do this without a loop?
    #    }
    color<-rep("black",n.points)
    color[n.points]<-"green"      # make the last one red
    cex<-rep(1.5,n.points) 
    cex[n.points]<-3
    hght.X<-history$height.X/(ifelse(max(history$height.X)<(10*pdf.max),10,(max(history$height.X)/pdf.max)))
    
    points(history$X,hght.X,pch="*",cex=cex,col=color)
  }
}


add.linear<-function(linear.all.out,num){
  for(i in 1:num){
    linear.all.out<-rbind(linear.all.out,do.one.linear(linear.all.out))    # need to loop this so that heights get added sequentially
  }
  return(linear.all.out)
  #  return(rbind(exp.all.out,rdply(num,do.one.exp(lambda,exp.all.out))[,-1]))  # [,-1] is to get rid of the first column (which is indices)
}

# most recent point details?

#################################################
# accept-reject

# function to run Accept-reject once
ar.runone<-function(fx,gx,Y,M,history){        
   U<-runif(1,0,1)
   fy<-fx(Y)
   gy<-gx(Y)
   status<-ifelse(U<=fy/(M*gy),"accept","reject")
   n.samples<-length(history$Y)
   if(n.samples>0){
     height.U<-sum(round(history$U[1:n.samples],digits=2)==round(U,digits=2))
     if(status=="accept"){
       height.Y<-sum(round(history$Y[1:n.samples],digits=2)==round(Y,digits=2) & history$status[1:n.samples]=="accept")  # can i do this without a loop?
     } else {
       height.Y<-(-sum(round(history$Y[1:n.samples],digits=2)==round(Y,digits=2) & history$status[1:n.samples]=="reject"))
     } 
   } else {
     height.U<-0
     height.Y<-0
   }
   return(data.frame(U=U,Y=Y,height.U,height.Y,status=status,fy=fy,gy=gy,M=M))
}

# Beta example
do.one.beta<-function(alpha,beta,all.out){
	gx<-function(x){return(dunif(x,0,1))}
	fx<-function(x){return(dbeta(x,shape1=alpha,shape2=beta))}
  y<-runif(1,0,1)    # use g(x)=Unif[0,1] for Beta example
	M<-optimize(dbeta,shape1=alpha,shape2=beta,interval=c(0,1),maximum=T)$objective
  # M is the maximum of the beta density function
	return(ar.runone(fx,gx,y,M,all.out))
}

plot.unif<-function(all.out){
  plot(NA,xlim=c(0,1),ylim=c(0,1),yaxt="n",xlab="u",ylab="")
  n.points<-length(all.out$U)
  if(n.points>0){
    cex<-rep(1.5,n.points)
    cex[n.points]<-3
    points(all.out$U,all.out$height.U/(ifelse(max(all.out$height.U)<10,10,max(all.out$height.U))),pch="*",cex=cex,col=c(rep("black",n.points-1),"green"))
  }
}

plot.beta<-function(alpha,beta,history){
  pdf.max<-optimize(dbeta,shape1=alpha,shape2=beta,interval=c(0,1),maximum=T)$objective
  x<-seq(0,1,by=0.01)
  plot(dbeta(x,alpha,beta)~x,type='l',xlab="y",ylab="density",ylim=c(0,pdf.max*1.5))
  if(alpha==1){
    lines(c(0,0),c(0,dbeta(0,alpha,beta)))
  }
  if(beta==1){
    lines(c(1,1),c(0,dbeta(1,alpha,beta)))
  }
  n.points<-dim(history)[1]
  if(n.points>0){
    hght.Y<-ifelse(history$status=="accept",history$height.Y/20,1.5*pdf.max+history$height.Y/20)
    if((max(history$height.Y)/20)>pdf.max){
     hght.Y<-ifelse(history$status=="accept",((pdf.max*history$height.Y)/(max(history$height.Y))),1.5*pdf.max+(1.5*pdf.max)*(history$height.Y/max(history$height.Y)))
    }
#    hght.Y<-history$height.Y/(ifelse(max(history$height.Y)<(10*pdf.max),10,(max(history$height.Y)/pdf.max)))
    color<-ifelse(history$status=="accept","black","gray80")
    color[n.points]<-ifelse(history$status[n.points]=="accept","green","red")      # make the last one red
    cex<-rep(1.5,n.points)
    cex[n.points]<-3
    points(history$Y,hght.Y,pch="*",cex=cex,col=color)
  }
}

temp.start<-function(alpha,beta){
  return(data.frame(U=NULL,V=NULL))
}

temp.start2<-function(beta.all.out,alpha,beta,num){
  for(i in 1:num){
    beta.all.out<-rbind(beta.all.out,do.one.beta(alpha,beta,beta.all.out))    # need to loop this so that heights get added sequentially
  }
  return(beta.all.out)
}


# Truncated Normal
do.one.tnorm<-function(all.out){
  gx<-function(x){return(exp(2-x))}   # might want to double-check this
  fx<-function(x){return(dnorm(x)/(1-pnorm(2)))}  # and this
  y<-2-log(1-runif(1,0,1))    # use PIT for f(y)=e^2e^-y1_{y \geq 2}
  M<-dnorm(2)/(1-pnorm(2))    
  return(ar.runone(fx,gx,y,M,all.out))
}

plot.tnorm<-function(history){
  pdf.max<-dnorm(2)/(1-pnorm(2))
  x<-seq(2,5,by=0.01)
  plot(dnorm(x)/(1-pnorm(2))~x,type='l',xlab="y",ylab="density",ylim=c(0,pdf.max*1.1))
  lines(pdf.max*exp(2-x)~x)

  n.points<-dim(history)[1]
  if(n.points>0){
    hght.Y<-ifelse(history$status=="accept",history$height.Y/20,pdf.max*exp(2-history$Y)+history$height.Y/20)
    if((max(history$height.Y)/20)>pdf.max){
      hght.Y<-ifelse(history$status=="accept",((pdf.max*history$height.Y)/(max(history$height.Y))),exp(2-history$Y)*pdf.max+(pdf.max)*(history$height.Y/max(history$height.Y)))
    }
    #    hght.Y<-history$height.Y/(ifelse(max(history$height.Y)<(10*pdf.max),10,(max(history$height.Y)/pdf.max)))
    color<-ifelse(history$status=="accept","black","gray80")
    color[n.points]<-ifelse(history$status[n.points]=="accept","green","red")      # make the last one red
    cex<-rep(1.5,n.points)
    cex[n.points]<-3
    points(history$Y,hght.Y,pch="*",cex=cex,col=color)
  }
}

start.tnorm<-function(){
  return(data.frame(NULL))
}

add.tnorm<-function(tnorm.all.out,num){
  for(i in 1:num){
    tnorm.all.out<-rbind(tnorm.all.out,do.one.tnorm(tnorm.all.out))    # need to loop this so that heights get added sequentially
  }
  return(tnorm.all.out)
}


#temp.start2<-function(linear.all.out,num){
#  for(i in 1:num){
##    linear.all.out<-rbind(linear.all.out,do.one.linear(linear.all.out))    # need to loop this so that heights get added sequentially
#  }
#  return(linear.all.out)
  #  return(rbind(exp.all.out,rdply(num,do.one.exp(lambda,exp.all.out))[,-1]))  # [,-1] is to get rid of the first column (which is indices)
#}


#add.one.beta<-function(alpha,beta){
#  
#}

server.R

library(shiny)
source("rvg.R")

shinyServer(function(input, output) {

#######################################
# Probability Integral Transform
#
 
# Exponential Example
  exp.all.out<-reactiveValues()
  observe({
    if(input$go.exp==0){exp.all.out$history<-init.all()} 
    else{
      exp.all.out$history<-add.exp(isolate(exp.all.out$history),isolate(input$lambda),isolate(input$num.exp))
    }
  })
  
  output$PITexpPlot <- renderPlot({
    input$go.exp
    input$clear.exp
    input$lambda
    par(mfrow=c(1,2),oma=c(0,0,0,0),mar=c(5.1,2.1,1,1.1))
    isolate(plot.unif(exp.all.out$history))
    isolate(plot.exp(input$lambda,exp.all.out$history))
  })
  
  observe({
    input$pitEx
    input$lambda
    input$clear.exp
    exp.all.out$history<-init.all()
  })
  

output$totalcountExp<-renderUI({
  input$go.exp
  last<-length(exp.all.out$history$X)
  if(last>0){
    isolate(paste("Total number of replicates: ",last))
  }
})


output$summaryExp <- renderUI({
  input$go.exp
  last<-length(exp.all.out$history$U)
  if(last>0){
    strexp1<-paste("The most recent value of U is:", 
                round(exp.all.out$history$U[last],3))
    strexp2<-"This gives the following for x:"
    HTML(paste(strexp1,strexp2,sep='<br/>'))
  }})


output$invExp<-renderUI({
  input$go.exp
  last<-length(exp.all.out$history$X)
  u<-exp.all.out$history$U[last]
  x<-exp.all.out$history$X[last]
  lambda<-input$lambda
  if(last>0){
      withMathJax(sprintf("$$x= \\frac{-ln(1-u)}{\\lambda} = \\frac{-ln(1-%0.3f)}{%0.1f} = %0.3f$$", u,lambda,x))
  }
})


# Linear example
linear.all.out<-reactiveValues()
observe({
  if(input$go.linear==0){linear.all.out$history<-init.all()} 
  else{
    linear.all.out$history<-add.linear(isolate(linear.all.out$history),isolate(input$num.linear))
  }
})

output$PITlinearPlot <- renderPlot({
  input$go.linear
  input$clear.linear
  par(mfrow=c(1,2),oma=c(0,0,0,0),mar=c(5.1,2.1,1,1.1))
  isolate(plot.unif(linear.all.out$history))
  isolate(plot.linear(linear.all.out$history))
})

observe({
  input$pitEx
  input$clear.linear
  linear.all.out$history<-init.all()
})


output$totalcountLin<-renderUI({
  input$go.linear
  last<-length(linear.all.out$history$X)
  if(last>0){
    isolate(paste("Total number of replicates: ",last))
  }
})



output$summaryLin <- renderUI({
  input$go.linear
  last<-length(linear.all.out$history$U)
  if(last>0){
    strexp1<-paste("The most recent value of U is:", 
                   round(linear.all.out$history$U[last],3))
    strexp2<-"This gives the following for x:"
    HTML(paste(strexp1,strexp2,sep='<br/>'))
  }})


output$invLin<-renderUI({
  input$go.linear
  last<-length(linear.all.out$history$X)
  u<-linear.all.out$history$U[last]
  x<-linear.all.out$history$X[last]
  if(last>0){
    withMathJax(sprintf("$$x= 4\\sqrt{u} = 4\\sqrt{%0.3f} = %0.3f$$", u,x))
  }
})

  
##########################################
# Accept-Reject 
#

# Beta example

  beta.all.out<-reactiveValues()
  observe({
    if(input$go==0){beta.all.out$history<-temp.start(isolate(input$alpha),isolate(input$beta))} 
    else{
    beta.all.out$history<-temp.start2(isolate(beta.all.out$history),isolate(input$alpha),isolate(input$beta),isolate(input$num))
  }
})

observe({
  input$alpha 
  input$beta 
  input$clear
  beta.all.out$history<-temp.start(input$alpha,input$beta)
  })


output$densityPlot <- renderPlot({
  input$go
  input$clear
  input$alpha
  input$beta
  par(mfrow=c(1,2),oma=c(0,0,0,0),mar=c(5.1,2.1,1,1.1))
  isolate(plot.unif(beta.all.out$history))
  isolate(plot.beta(input$alpha,input$beta,beta.all.out$history))
})

output$summary <- renderUI({
    input$go
    last<-length(beta.all.out$history$Y)
    if(last>0){
      str1<-paste("The most recent value of U is:", 
                  round(beta.all.out$history$U[last],3), "(enlarged and green)")
      str2<-paste("The most recent value of Y is:", 
            round(beta.all.out$history$Y[last],3), "(enlarged, and green if accepted; red if rejected)")
      str3<-paste("The value of Y is", ifelse(beta.all.out$history$status[last]=="accept","<b>accepted</b>","<b>rejected</b>"),"because:")
      HTML(paste(str1,str2,str3,sep='<br/>'))
}})

output$accrej<-renderUI({
  input$go
  last<-length(beta.all.out$history$Y)
  if(last>0){
    u<-beta.all.out$history$U[last]
    fy<-beta.all.out$history$fy[last]
    M<-beta.all.out$history$M[last]
    gy<-beta.all.out$history$gy[last]
    ratio<-fy/(M*gy)
    if(beta.all.out$history$status[last]=="accept"){
       withMathJax(sprintf("$$%0.3f \\leq \\frac{f(y)}{Mg(y)} = 
\\frac{\\frac{\\Gamma(\\alpha + \\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)}y^{\\alpha-1}(1-y)^{\\beta-1}}{M \\cdot 1_{\\{0 \\leq y \\leq 1\\}}}
                           = \\frac{%0.2f}{%0.3f \\cdot 1} = %0.2f$$",
                           u,fy,M,ratio))
     } else {
       withMathJax(sprintf("$$%0.3f > \\frac{f(y)}{Mg(y)} = 
\\frac{\\frac{\\Gamma(\\alpha + \\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)}y^{\\alpha-1}(1-y)^{\\beta-1}}{M \\cdot 1_{\\{0 \\leq y \\leq 1\\}}}
                           = \\frac{%0.2f}{%0.2f \\cdot 1} = %0.2f$$",
                           u,fy,M,ratio))
     }
  }
})

output$unifnote<-renderText({
  input$alpha
  input$beta
  if(input$alpha==1 & input$beta==1){
    "Note that for the Beta(1,1) distribution, every point will be accepted, as we would expect
    since it is equivalent to the Uniform[0,1] distribution."
  }
  
})

output$M<- renderText({
  input$alpha
  input$beta
  isolate(paste("For the current set of parameter values, M = ", 
                  round(optimize(dbeta,shape1=input$alpha,shape2=input$beta,interval=c(0,1),maximum=T)$objective,digits=3),".",sep=""))
})

output$totalcount<-renderUI({
  input$go
  last<-length(beta.all.out$history$Y)
  if(last>0){
    isolate(paste("Total number of replicates: ",last))
  }
})


# Truncated normal example

tnorm.all.out<-reactiveValues()
observe({
  if(input$tnormGo==0){tnorm.all.out$history<-start.tnorm()} 
  else{
    tnorm.all.out$history<-add.tnorm(isolate(tnorm.all.out$history),isolate(input$tnormNum))
  }
})

observe({
  input$tnormClear
  tnorm.all.out$history<-start.tnorm()
})



output$tnormDensityPlot <- renderPlot({
  input$tnormGo
  input$tnormClear
  par(mfrow=c(1,2),oma=c(0,0,0,0),mar=c(5.1,2.1,1,1.1))
  isolate(plot.unif(tnorm.all.out$history))
  isolate(plot.tnorm(tnorm.all.out$history))
})

output$tnormsummary <- renderUI({
  input$tnormGo
  last<-length(tnorm.all.out$history$Y)
  if(last>0){
    str1<-paste("The most recent value of U is:", 
                round(tnorm.all.out$history$U[last],3), "(enlarged and green)")
    str2<-paste("The most recent value of Y is:", 
                round(tnorm.all.out$history$Y[last],3), "(enlarged, and green if accepted; red if rejected)")
    str3<-paste("The value of Y is", ifelse(tnorm.all.out$history$status[last]=="accept","<b>accepted</b>","<b>rejected</b>"),"because:")
    HTML(paste(str1,str2,str3,sep='<br/>'))
  }})


output$tnormaccrej<-renderUI({
  input$tnormGo
  last<-length(tnorm.all.out$history$Y)
  if(last>0){
    u<-tnorm.all.out$history$U[last]
    fy<-tnorm.all.out$history$fy[last]
    M<-tnorm.all.out$history$M[last]
    gy<-tnorm.all.out$history$gy[last]
    y<-tnorm.all.out$history$Y[last]
    ratio<-fy/(M*gy)
    if(tnorm.all.out$history$status[last]=="accept"){
      withMathJax(sprintf("$$%0.3f \\leq \\frac{f(y)}{Mg(y)} = 
                          \\frac{\\frac{1}{\\sqrt{2 \\pi}}e^{-\\frac{1}{2}(%0.2f)^2}
                    \\cdot \\left[\\frac{1}{1-\\Phi(2)}\\right] }{M \\cdot e^{2-%0.2f} }
                          = \\frac{%0.2f}{%0.3f \\cdot %0.3f} = %0.2f$$",
                          u,y,y,fy,M,gy,ratio))
    } else {
      withMathJax(sprintf("$$%0.3f > \\frac{f(y)}{Mg(y)} = 
                          \\frac{\\frac{1}{\\sqrt{2 \\pi}}e^{-\\frac{1}{2}(%0.2f)^2}
                    \\cdot \\left[\\frac{1}{1-\\Phi(2)}\\right] }{M \\cdot e^{2-%0.2f} }
                          = \\frac{%0.2f}{%0.3f \\cdot %0.3f} = %0.2f$$",
                          u,y,y,fy,M,gy,ratio))
    }
  }
  })

output$tnormRatio<-renderUI({
  input$tnormGo
  num<-length(tnorm.all.out$history$Y)
  if(num>0){
    str4<-paste("The proportion of points that have been accepted is <b>", 
                round(sum(tnorm.all.out$history$status=="accept")/num,3),"</b> (out of ",num,")",sep="")
    HTML(str4)    
  }
})



})  # end of shinyServer

ui.R

# ---------------------------------------
#  App Title: Random Variable Generation
#     Author: Peter Chi
# ---------------------------------------

library(shiny)
source("rvg.R")

shinyUI(navbarPage("Random Variable Generation",
    tabPanel("(1) Probability Integral Transform",
    fluidPage(withMathJax(),
	       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")),  
      h3("Random Variable Generation: Probability Integral Transform"),
      p("Suppose we would like to generate \\(X\\sim f\\), where \\(f\\) is the probability density function (pdf) of \\(X\\). If the corresponding
        cumulative distribution function (cdf) has a generalized inverse, then we can use the Probability Integral Transform. The only other requirement
        is that we have the ability to simulate \\(U\\sim Unif[0,1]\\)."),
          sidebarLayout(
       sidebarPanel(
         selectizeInput('pitEx', label=("Choose an example:"), choices=c(Exponential="Exponential",Other="Other")),
         conditionalPanel(
           condition="input.pitEx == 'Exponential'",
           sliderInput(inputId="lambda",label="Value of \\(\\lambda\\) (rate) parameter:",min=0.1,max=5,value=1,step=0.1),
           br(),br(),
           sliderInput(inputId="num.exp",label="Generate how many?",min=1,max=500,value=1,step=1),
           div(actionButton("clear.exp",label="Clear"),actionButton("go.exp",label="Generate"),align="right"),
           br(),
           uiOutput("totalcountExp")
         ),
         conditionalPanel(
           condition="input.pitEx == 'Other'",
           helpText("This example is with an arbitrary, un-named distribution (i.e.
                    one for which pre-packaged routines are unlikely to exist)."),
           br(),br(),
           sliderInput(inputId="num.linear",label="Generate how many?",min=1,max=500,value=1,step=1),
           div(actionButton("clear.linear",label="Clear"),actionButton("go.linear",label="Generate"),align="right"),
           br(),
           uiOutput("totalcountLin")
         ),
         p(),br(),br(),
         
         div("Shiny app by", a(href= "http://statweb.calpoly.edu/pchi/", target="_blank", "Peter Chi"),align="right", style = "font-size: 8pt"),
         div("Base R code by", a(href= "http://statweb.calpoly.edu/pchi/", target="_blank", "Peter Chi"),align="right", style = "font-size: 8pt"),
         div("Shiny source files:", a(href= "https://gist.github.com/calpolystat/cc55e764b757d8729963", 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")
       ),
       mainPanel(
         conditionalPanel(
                condition = "input.pitEx == 'Exponential'",  
                h3("Demonstration with \\(X \\sim Exp(\\lambda)\\):"),
                plotOutput("PITexpPlot",width="100%"),
 #               helpText("Most recent point is enlarged and green."),
 uiOutput("summaryExp"),
 uiOutput("invExp"),
                HTML("<hr style='height: 2px; color: #de7008; background-color: #df7109; border: none;'>"),
                p("Details:"),
                p("1) First we note that as \\(f(x)=\\lambda e^{-\\lambda x}\\) for an Exponential random variable, 
                  the cdf is thus \\(F(x) = 1-e^{-\\lambda x}\\), for \\(x \\geq 0\\)"),
                p("2) Next, the inverse of this is \\(F^{-1}(u) = \\frac{-ln(1-u)}{\\lambda}\\)"),
                p("3) Thus, we generate \\(U\\sim Unif[0,1]\\), and plug these values into \\(F^{-1}\\) to obtain generations of \\(X\\)")
          ),
          
          conditionalPanel(
            condition = "input.pitEx == 'Other'", 
            h3("Demonstration with \\(f(x) = \\frac{x}{8} \\cdot 1_{\\{0 \\leq x \\leq 4\\}}\\) "),
            plotOutput("PITlinearPlot",width="100%"),
#            helpText("Most recent point is enlarged and green."),
uiOutput("summaryLin"),
uiOutput("invLin"),
            HTML("<hr style='height: 2px; color: #de7008; background-color: #df7109; border: none;'>"),
            p("Details:"),
            p("1) First we note that the cdf is \\(F(x) = \\frac{x^2}{16}\\), for \\(0 \\leq x \\leq 4\\)"),
            p("2) Next, the inverse of this is \\(F^{-1}(u) = 4 \\sqrt{u}\\)"),
            p("3) Thus, we generate \\(U\\sim Unif[0,1]\\), and plug these values into \\(F^{-1}\\) to obtain generations of \\(X\\)")
        #    plotOutput("PITexpPlot",width="100%")
          )
#          )
)
)
)
),

    tabPanel("(2) Accept-Reject Algorithm",
             fluidPage(withMathJax(),
                       h3("Random Variable Generation: Accept-Reject Algorithm"),
                      p("Suppose we would like to generate \\(X\\sim f\\) but can neither do it directly
nor via the Probability Integral Transform (e.g. if the generalized inverse
of the cdf is unavailable). We can instead arrive at it via generating \\(Y\\sim g\\),
with only the following two necessary conditions:"),
                       p("1) \\(f\\) and \\(g\\) have the same support"),
                      p("2) We can find a constant \\(M\\) such that \\(f(x)/g(x) \\leq M \\) for all \\(x\\) "),
               sidebarLayout(
                 sidebarPanel(
                   selectInput("example", label=("Choose an example:"), choices=list("Beta"="Beta","Truncated Normal"="tnorm")),
                 conditionalPanel(
                   condition = "input.example == 'Beta'",
                  sliderInput(inputId="alpha",label="Value of \\(\\alpha\\) parameter:",min=1,max=5,value=1,step=0.1),
                   sliderInput(inputId="beta",label="Value of \\(\\beta\\) parameter:",min=1,max=5,value=1,step=0.1),
                   br(),br(),
                   sliderInput(inputId="num",label="Generate how many?",min=1,max=500,value=1,step=1),
                   div(actionButton("clear",label="Clear"),actionButton("go",label="Generate"),align="right"),
                  br(),
                  uiOutput("totalcount")
                 ),
                 conditionalPanel(
                   condition = "input.example == 'tnorm'",
                   br(),
                   sliderInput(inputId="tnormNum",label="Generate how many?",min=1,max=500,value=1,step=1),
                   div(actionButton("tnormClear",label="Clear"),actionButton("tnormGo",label="Generate"),align="right"),
                   br(),br(),
                   helpText("This example is with the standard normal distribution, truncated at 2 (i.e. allowing for values greater than or equal to 2 only)."),
                   br(),
                   helpText("It would indeed be possible to simulate a normal random variable truncated at 2 by using a pre-packaged
                            routine for a standard normal random variable (such as rnorm in R), and then discarding all values below 2."),
                  br(),
                  helpText("However, this would be extremely inefficient as we would discard more than 97% of all values that we generate.
                            With the accept-reject algorithm, we of course also discard values, but here we demonstrate that this method
                            has superior efficiency, in the sense that less than 97% of the generated values will be discarded.")
                   ),
                 
                 p(),br(),br(),
                 
                 div("Shiny app by", a(href= "http://statweb.calpoly.edu/pchi/", target="_blank", "Peter Chi"),align="right", style = "font-size: 8pt"),
                 div("Base R code by", a(href= "http://statweb.calpoly.edu/pchi/", target="_blank", "Peter Chi"),align="right", style = "font-size: 8pt"),
                 div("Shiny source files:", a(href= "https://gist.github.com/calpolystat/cc55e764b757d8729963", 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(7,
                 conditionalPanel(
                   condition = "input.example == 'Beta'",                   
                   h3("Demonstration with \\(X \\sim Beta(\\alpha,\\beta)\\):"),
                   plotOutput("densityPlot",width="100%"),
                   uiOutput("summary"),
                   uiOutput("accrej"),
                   br(),
                   helpText("In the right panel: after initially being shown in red, rejected points remain in grey and stack down from the top; after initially
being shown in green, accepted
                            points remain in black and stack up from the bottom, to fill the shape of the theoretical pdf of \\(X\\)."),
                   textOutput("unifnote"),
                   br(),
                   HTML("<hr style='height: 2px; color: #de7008; background-color: #df7109; border: none;'>"),
                   p("Details:"),
                   p("1) The first step is to generate \\(Y \\sim g\\). In this example, we use \\(Y \\sim Unif[0,1]\\) 
                     (shown in the right panel, along with the true distribution that we are trying to simulate from). Notice
                     that the Unif[0,1] distribution does indeed have the same support as the Beta distribution."),
                   p("2) Next, we need to find an appropriate value of M. For the Beta example, 
                     we notice that the maximum of the Beta pdf would work."),
                   textOutput("M"),
                   br(),
                   p("3) We also generate \\(U \\sim Unif[0,1]\\) (left panel), and then accept \\(y\\) as a 
                     value of \\(X\\) if \\(U \\leq \\frac{f(y)}{Mg(y)}\\), and reject otherwise")
                 ),
                 conditionalPanel(
                   condition = "input.example == 'tnorm'",
                   h3("Demonstration with Truncated Normal"),
                   plotOutput("tnormDensityPlot",width="100%"),
                   uiOutput("tnormsummary"),
                   uiOutput("tnormaccrej"),
                   br(),
                   helpText("In the right panel: after initially being shown in red, rejected points remain in grey and stack down from the top; after initially
being shown in green, accepted
                            points remain in black and stack up from the bottom, to fill the shape of the theoretical pdf of \\(X\\)."),
                   br(),
                   uiOutput("tnormRatio"),
                   br(),
                 HTML("<hr style='height: 2px; color: #de7008; background-color: #df7109; border: none;'>"),
                 p("Details:"),
                 p("1) The first step is to generate \\(Y \\sim g\\). In this example, we will use \\(g(y) = e^{2-y} \\cdot 1_{\\{y \\geq 2\\}} \\) 
                     (shown in the right panel, along with the true distribution that we are trying to simulate from)."),
                 p("2) Next, we need to find an appropriate value of M. For this example, 
                     we notice the following:"),
                 p("$$\\frac{f(y)}{g(y)} = \\frac{\\frac{1}{\\sqrt{2 \\pi}}e^{-\\frac{1}{2}y^2}1_{\\{y \\geq 2\\}}
                    \\cdot \\left[\\frac{1}{1-\\Phi(2)}\\right] }{e^{2-y}1_{\\{y \\geq 2\\}}}$$"),
                 p("where \\(\\Phi\\) is the standard normal cdf.
                   It can be shown that this ratio is at its maximum at \\(y=2\\). Thus, \\(M=\\frac{\\phi(2)}{1-\\Phi(2)}\\) where \\(\\phi\\) is
                   the standard normal pdf."),
                 p("3) We also generate \\(U \\sim Unif[0,1]\\) (left panel), and then accept \\(y\\) as a 
                     value of \\(X\\) if \\(U \\leq \\frac{f(y)}{Mg(y)}\\), and reject otherwise")                 
                 )
                 )
               )
             )
    )
)
)