Last active
August 29, 2015 14:23
-
-
Save calpolystat/cc55e764b757d8729963 to your computer and use it in GitHub Desktop.
Revisions
-
calpolystat revised this gist
Jun 20, 2015 . 2 changed files with 2 additions and 3 deletions.There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -1,4 +1,3 @@ ############################################# # Probability Integral Transform # This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -40,7 +40,7 @@ shinyUI(navbarPage("Random Variable Generation", 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( @@ -121,7 +121,7 @@ with only the following two necessary conditions:"), 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") -
Jun 20, 2015 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,7 @@ 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 This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,7 @@ 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 This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,21 @@ 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. This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,245 @@ ############################################# # 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){ # #} This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,294 @@ 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 This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,186 @@ # --------------------------------------- # 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/96f9adc5c37f4414fbd1", 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/96f9adc5c37f4414fbd1", 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") ) ) ) ) ) ) )