R sampling app version 3

Continuing from my previous post, R sampling app version 2 [app], here I expand upon the example Shiny app that was presented, which generates random samples from some common probability distributions. Specifically, I incorporate additional probability density and mass functions.

RshinyIntro3
You can click the image to go to the app page.

See my posts, Mathematical notation in R plots and Mathematical notation in R plots 2, for some examples and conveniently supplied R code for displaying the equation of a pdf or pmf on an R plot for many common distributions. The code shown in this post for version 3 of the random sampling app will make use of all these R expressions, but you won’t see the expressions explicitly, so please see the other posts if you are interested in how they work or just copying them for your convenience.

I have done some cleanup of the scripts. First, here is the code that is run to build the workspace file that will be loaded by the app:

n.default <- 500
#discrete r* wrapper functions
rbern <- function(n=n.default,bern.prob=0.5){ rbinom(n=n,size=1,prob=bern.prob) }
rbinom2 <- function(n=n.default,binom.size=10,binom.prob=0.5){ rbinom(n,size=binom.size,prob=binom.prob) }
drunif <- function(n=n.default,drunif.min=0,drunif.max=100,drunif.step=1){ sample(seq(drunif.min,drunif.max,by=drunif.step),size=n,rep=T) }
rgeom2 <- function(n=n.default,geom.prob=0.5){ rgeom(n,prob=geom.prob) }
rhyper2 <- function(n=n.default,hyper.M=10,hyper.N=20,hyper.K=10){ rhyper(nn=n,m=hyper.M,n=hyper.N-hyper.M,k=hyper.K) }
rnbinom2 <- function(n=n.default,nbin.size=10,nbin.prob=0.5){ rnbinom(n,size=nbin.size,prob=nbin.prob) }
rpois2 <- function(n=n.default,poi.lambda=10){ rpois(n,poi.lambda) }
# continuous r* wrapper functions
rbeta2 <- function(n=n.default,beta.shape1=2,beta.shape2=2){ rbeta(n,shape1=beta.shape1,shape2=beta.shape2) }
rcauchy2 <- function(n=n.default,cau.location=0,cau.scale=1){ rcauchy(n,location=cau.location,scale=cau.scale) }
rchisq2 <- function(n=n.default,chisq.df=1){ rchisq(n,df=chisq.df) }
rexp2 <- function(n=n.default,exp.rate=1){ rexp(n=n,rate=exp.rate) }
rf2 <- function(n=n.default,F.df1=1,F.df2=15){ rf(n,df1=F.df1,df2=F.df2) }
rgamma2 <- function(n=n.default,gam.shape=1,gam.rate=1){ rgamma(n,shape=gam.shape,rate=gam.rate) }
rlaplace2 <- function(n=n.default,lap.location=0,lap.scale=1){ rlaplace(n,location=lap.location,scale=lap.scale) }
rlogis2 <- function(n=n.default,logi.location=0,logi.scale=1){ rlogis(n,location=logi.location,scale=logi.scale) }
rpareto2 <- function(n=n.default,pareto.location=1,pareto.shape=3){ rpareto(n,location=pareto.location,shape=pareto.shape) }
rweibull2 <- function(n=n.default,weib.shape=1,weib.scale=1){ rweibull(n,shape=weib.shape,scale=weib.scale) }
rt2 <- function(n=n.default,t.df=15){ rt(n=n,df=t.df) }

# Continuous distribution plotmath expressions:
expr.beta <- expression(italic(paste(displaystyle(f(x)~"="~frac(Gamma(alpha+beta),Gamma(alpha)*Gamma(beta))*x^{alpha-1}*(1-x)^{beta-1})
					~~~~displaystyle(list(paste(0<=x) <=1, atop(paste(0<alpha) <infinity, paste(0<beta) <infinity)))
					)))

expr.cauchy <- expression(italic(paste(displaystyle(f(x)~"="~frac(1,pi*sigma)~frac(1,1+bgroup("(",frac(x-theta,sigma),")")^2))
					~~~~displaystyle(list(paste(-infinity<x) <infinity, atop(paste(-infinity<theta) <infinity, sigma > 0)))
					)))

expr.chisq <- expression(italic(paste(frac(1,2^{frac(nu,2)}*Gamma~bgroup("(",frac(nu,2),")"))*x^{frac(nu,2)-1}*e^{-frac(x,2)}
					~~~~displaystyle(atop(paste(0<=x) <infinity, nu~"="~list(1,2,...)))
					)))

expr.exp <- expression(italic(paste(displaystyle(f(x)~"="~lambda*e^{-lambda*x})
					~~~~displaystyle(atop(paste(0<=x) <infinity,lambda>0))
					)))

expr.F <- expression(italic(paste(displaystyle(f(x)~"="~frac(Gamma~bgroup("(",frac(nu[1]+nu[2],2),")"),Gamma~bgroup("(",frac(nu[1],2),")")~Gamma~bgroup("(",frac(nu[2],2),")"))
					~bgroup("(",frac(nu[1],nu[2]),")")^{frac(nu[1],2)}~frac(x^{frac(nu[1],2)-1},bgroup("(",1+frac(d[1],d[2])*x,")")^{frac(d[1]+d[2],2)}))
					~~~~displaystyle(atop(paste(0<=x) <infinity,list(d[1],d[2])~"="~list(1,2,...)))
					)))

expr.gam <- expression(italic(paste(displaystyle(f(x)~"="~frac(beta^alpha,Gamma(alpha))*x^{alpha-1}*e^{-beta*x})
					~~~~displaystyle(list(paste(0<x) <infinity, atop(paste(0<alpha) <infinity, paste(0<beta) <infinity)))
					)))

expr.lap <- expression(italic(paste(displaystyle(f(x)~"="~frac(1,2*sigma)~e^{-frac(abs(x-mu),sigma)})
					~~~~displaystyle(list(paste(-infinity<x) <infinity, atop(paste(-infinity<mu) <infinity, sigma > 0)))
					)))

expr.logi <- expression(italic(paste(displaystyle(f(x)~"="~frac(1,beta)~frac(e^{-frac(x-mu,beta)},bgroup("[",1+e^{-frac(x-mu,beta)},"]")^2))
					~~~~displaystyle(list(paste(-infinity<x) <infinity, atop(paste(-infinity<mu) <infinity, beta > 0)))
					)))

expr.lognorm <- expression(italic(paste(displaystyle(f(x)~"="~frac(1,x*sigma*sqrt(2*pi))~e^{-frac((log(x)-mu)^2,2*sigma^2)})
					~~~~displaystyle(list(paste(0<x) <infinity, atop(paste(-infinity<log(mu)) <infinity, paste(0<sigma^scriptscriptstyle("2")) <infinity)))
					)))

expr.norm <- expression(italic(paste(displaystyle(f(x)~"="~frac(1,sqrt(2*pi*sigma^scriptscriptstyle("2")))~e^{frac(-1,2*sigma^{scriptscriptstyle("2")})*(x-mu)^scriptscriptstyle("2")})
					~~~~displaystyle(list(paste(-infinity<x) <infinity, atop(paste(-infinity<mu) <infinity, paste(0<sigma^scriptscriptstyle("2")) <infinity)))
					)))

expr.pareto <- expression(italic(paste(displaystyle(f(x)~"="~frac(beta*alpha^beta,x^{beta+1}))
					~~~~displaystyle(atop(paste(alpha<x) <infinity, list(alpha,beta) > 0))
					)))

expr.t <- expression(italic(paste(displaystyle(f(x)~"="~frac(Gamma~bgroup("(",frac(nu+1,2),")"),sqrt(nu*pi)~Gamma~bgroup("(",frac(nu,2),")"))~bgroup("(",1+frac(x^2,nu),")")^{-frac(nu+1,2)})
					~~~~displaystyle(atop(paste(-infinity<x) <infinity, nu > 0))
					)))

expr.unif <- expression(italic(paste(displaystyle(f(x)~"="~frac(1,b-a)
					~~~~displaystyle(paste(-infinity<paste(a<=paste(x<=b))) <infinity))
					)))

expr.weib <- expression(italic(paste(displaystyle(f(x)~"="~frac(k,lambda)~bgroup("(",frac(x,lambda),")")^{k-1}*e^(-x/lambda)^k)
					~~~~displaystyle(atop(paste(0<=x) <infinity, list(k,lambda) > 0))
					)))

# Discrete distribution plotmath expressions:
expr.bern <- expression(italic(paste(displaystyle(P(X~"="~x)~"="~p^x*(1-p)^{1-x})
					~~~~displaystyle(atop(x~"="~list(0,1), paste(0<=p)<=1))
					)))

expr.bin <- expression(italic(paste(displaystyle(P(X~"="~x)~"="~bgroup("(",atop(n,x),")")~p^x*(1-p)^{n-x})
					~~~~displaystyle(atop(x~"="~list(0,1,...,n), paste(0<=p)<=1))
					)))

expr.dunif <- expression(italic(paste(displaystyle(P(X~"="~x)~"="~frac(1,N))
					~~~~displaystyle(x~"="~list(1,2,...,N))
					)))

expr.geom <- expression(italic(paste(displaystyle(P(X~"="~x)~"="~p*(1-p)^x)
					~~~~displaystyle(atop(x~"="~list(1,2,...), paste(0<=p)<=1))
					)))

expr.hgeom <- expression(italic(paste(displaystyle(P(X~"="~x)~"="~frac(bgroup("(",atop(M,x),")")~bgroup("(",atop(N-M,K-x),")"),bgroup("(",atop(N,K),")")))
					~~~~displaystyle(list(x~"="~list(0,1,...,K), atop(paste(M-(N-K)<=x)<=M, list(N,M,K)>=0)))
					)))

expr.nbin <- expression(italic(paste(displaystyle(P(X~"="~x)~"="~frac(Gamma(x+n),Gamma(n)*x*"!")~p^r*(1-p)^x)
					~~~~displaystyle(atop(x~"="~list(0,1,...), paste(0<=p)<=1))
					)))

expr.poi <- expression(italic(paste(displaystyle(P(X~"="~x)~"="~frac(e^{-lambda}*lambda^x,x*"!"))
					~~~~displaystyle(atop(x~"="~list(0,1,...), paste(0<=lambda)<infinity))
					)))

I save this workspace as samplingApp.RData and place it in the directory that contains the ui.R and server.R scripts. Next, we have the ui.R script. The most important, though subtle, change is the addition of uiOutput("dist3") on account of now including some distributions in the app which take three input arguments when sampling from them. You can see here all 21 of the distributions, 7 discrete and 14 continuous, which have been incorporated into the app.

shinyUI(pageWithSidebar(
	headerPanel("Distributions of Random Variables"),
	sidebarPanel(
		radioButtons("dist","Distribution type:",
			list(
				"Bernoulli"="bern","Binomial"="bin","Discrete Uniform"="dunif","Geometric"="geom","Hypergeometric"="hgeom","Negative Binomial"="nbin","Poisson"="poi", # discrete
				"Beta"="beta","Cauchy"="cauchy","Chi-squared"="chisq","Exponential"="exp","F"="F","Gamma"="gam","Laplace (Double Exponential)"="lap", # continuous
				"Logistic"="logi","Log-Normal"="lnorm","Normal"="norm","Pareto"="pareto","t"="t","Uniform"="unif","Weibull"="weib"
				)
		),
		sliderInput("n","Sample size:",1,1000,500),
		uiOutput("dist1"),
		uiOutput("dist2"),
		uiOutput("dist3"),
		checkboxInput("density","Show density curve",FALSE),
		conditionalPanel(
			condition="input.density==true",
			numericInput("bw","bandwidth:",1)
		),
		downloadButton('dldat', 'Download Sample')
	),
	mainPanel(
		tabsetPanel(
			tabPanel("Plot",plotOutput("plot",height="auto")),
			tabPanel("Summary",verbatimTextOutput("summary")),
			tabPanel("Table",tableOutput("table"))
		)
	)
))

Finally we have the server.R script, which has changed significantly since version 2. I have cleaned it up by alphabetizing things as the number of distributions has grown. Also, I moved much of the code in the top of the script prior to the shinyServer call to the workspace file mentioned above. You may have noted there, if you saw my previous posts, that I removed any references to the formals command and standardized all of my r* random sampling wrapper functions. The reason I made these wrapper functions in the first place, e.g. rgamma2 instead of just using the original base R functions, was because the reactivity of Shiny made it a challenge to work with multiple functions that had formal arguments of the same name. Perhaps there is a better way around this, but this was what I thought to do. Corresponding to the change noted for the ui.R script, in the server.R script I have added output$dist3. Lastly, the VGAM package is loaded to allow sampling from the Laplace and Pareto distributions.

library(shiny)
pkgs <- c("VGAM")
pkgs <- pkgs[!(pkgs %in% installed.packages()[,"Package"])]
if(length(pkgs)) install.packages(pkgs,repos="http://cran.cs.wwu.edu/")
library(VGAM)
load("samplingApp.RData", envir=.GlobalEnv)

shinyServer(function(input,output){
	dat <- reactive({
		dist <- switch(input$dist,
			bern=rbern, bin=rbinom2, dunif=drunif, geom=rgeom2, hgeom=rhyper2, nbin=rnbinom2, poi=rpois2, # discrete
			beta=rbeta2, cauchy=rcauchy2, chisq=rchisq2, exp=rexp2, F=rf2, gam=rgamma2, lap=rlaplace2, # continuous
			logi=rlogis2, lnorm=rlnorm, norm=rnorm, pareto=rpareto2, t=rt2, unif=runif, weib=rweibull2,
			)

		def.args <- switch(input$dist,
			# discrete
			bern=c(input$bern.prob),
			bin=c(input$binom.size,input$binom.prob),
			dunif=c(input$drunif.min,input$drunif.max,input$drunif.step),
			geom=c(input$geom.prob),
			hgeom=c(input$hyper.M,input$hyper.N,input$hyper.K),
			nbin=c(input$nbin.size,input$nbin.prob),
			poi=c(input$poi.lambda),
			# continuous
			beta=c(input$beta.shape1,input$beta.shape2),
			cauchy=c(input$cau.location,input$cau.scale),
			chisq=c(input$chisq.df),
			exp=c(input$exp.rate),
			F=c(input$F.df1,input$F.df2),
			gam=c(input$gam.shape,input$gam.rate),
			lap=c(input$lap.location,input$lap.scale),
			logi=c(input$logi.location,input$logi.scale),
			lnorm=c(input$meanlog,input$sdlog),
			norm=c(input$mean,input$sd),
			pareto=c(input$pareto.location,input$pareto.shape),
			t=c(input$t.df),
			unif=c(input$min,input$max),
			weib=c(input$weib.shape,input$weib.scale)
			)

		f <- formals(dist)
		f <- f[names(f)!="nn" & names(f)!="n"]
		if(any(input$dist==c("dunif","hgeom"))){ len <- min(length(f),4-1); f <- f[1:len]
		} else { len <- min(length(f),3-1); f <- f[1:len] }
		argList <- list(n=input$n)
		for(i in 1:len) argList[[names(f)[i]]] <- def.args[i]
		return(list(do.call(dist,argList),names(f)))
	})

	output$dist1 <- renderUI({
		lab <- switch(input$dist,
			bern="Probability:", bin="Size:", dunif="Discrete sequence minimum:", geom="Probability:", hgeom="M:", nbin="Number of successes:",	poi="Mean and Variance:", # discrete
			beta="Alpha:", cauchy="Location:", chisq="Degrees of freedom:", exp="Rate:", F="Numerator degrees of freedom:", gam="Shape:", lap="Location:",
			logi="Location:", lnorm="Mean(log):", norm="Mean:", pareto="Location:",	t="Degrees of freedom:", unif="Minimum:", weib="Shape:"
			# logistic and log-normal are switched on plotmath expressions graphic!
			)
		ini <- switch(input$dist,
			bern=0.5, bin=10, dunif=0, geom=0.5, hgeom=10, nbin=10, poi=10,	# discrete
			beta=2, cauchy=0, chisq=1, exp=1, F=1, gam=1, lap=0, logi=0, lnorm=0, norm=0, pareto=1,	t=15, unif=0, weib=1 # continuous
			)
		numericInput(dat()[[2]][1],lab,ini)
	})
	
	output$dist2 <- renderUI({
		lab <- switch(input$dist,
			bin="Probability:",	dunif="Discrete sequence maximum:",	hgeom="N:",	nbin="Probability:", # discrete
			beta="Beta:", cauchy="Scale:", F="Denominator degrees of freedom:", gam="Rate:", lap="Scale:", # continuous
			logi="Scale:", lnorm="Standard deviation(log)", norm="Standard deviation:", pareto="Shape:", unif="Maximum:", weib="Scale:"
			)
		ini <- switch(input$dist,
			bin=0.5, dunif=100, hgeom=20, nbin=0.5,
			beta=2, cauchy=1, F=15, gam=1, lap=1, logi=1, lnorm=1, norm=1, pareto=3, unif=1, weib=1
			)
		if(any(input$dist==c("bin","dunif","hgeom","nbin","cauchy","lap","logi","pareto","weib",
							"beta","F","gam","lnorm","norm","unif"))) numericInput(dat()[[2]][2],lab,ini)
	})
	
	output$dist3 <- renderUI({
		lab <- switch(input$dist,
			dunif="Step size:",	hgeom="K:")
		ini <- switch(input$dist,
			dunif=1, hgeom=5)
		if(any(input$dist==c("dunif","hgeom"))) numericInput(dat()[[2]][3],lab,ini)
	})
	
	output$dldat <- downloadHandler(
		filename = function() { paste(input$dist, '.csv', sep='') },
		content = function(file) {
			write.csv(data.frame(x=dat()[[1]]), file)
		}
	)

	output$plot <- renderPlot({
		dist <- input$dist
		n <- input$n
		expr <- get(paste("expr",dist,sep="."))
		par(mar=c(4,4,10,1))
		hist(dat()[[1]],main=expr,xlab="Observations",col="orange",cex.main=1.5,cex.axis=1.3,cex.lab=1.3,prob=T)
		if(input$density) lines(density(dat()[[1]],adjust=input$bw),lwd=2)
	},
	height=750, width=1000
	)
	
	output$summary <- renderPrint({
		summary(dat()[[1]])
	})
	
	output$table <- renderTable({
		data.frame(x=dat()[[1]])
	})
	
})

Go to version 3 of the sampling app. Here is the blog post for version 4.

This entry was posted by Matt Leonawicz.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: