#' Run the Probability Walkthrough script.
#'
#' @export
#' @examples
#' \dontrun{prob()}
prob<-function()
{
##' @author Josh Errickson
getvalidinput <- function(prompt, error, validinput = NULL, low=0, high=0)
{
if(is.null(validinput) && low == high) stop("Must either input validinput or both low and high!")
while(TRUE)
{
input <- readline(prompt)
if(input %in% validinput)
{
return(input)
}
if(!is.na(as.numeric(input)))
{
if(low < high && low < high && low <= as.numeric(input) && as.numeric(input) <= high)
{
return(input)
}
}
cat(error)
}
}
options(warn=-1) # make this 0 for testing, -1 for production
options(scipen=4) #Don't give numbers like 1e-4, instead do .0001
cat("Are you ready for some nice graphics to help illustrate probabilities?\n")
cat("At each step, type your response and press \"Enter\"\n")
while(TRUE)
{
cat("What is the form of your probability? Choose from the following:\n1: P(Z < c)\n2: P(Z > c)\n3: P(a < Z < b)\n4: P(Z < a or Z > b)\nYou can press 'q' at any input to quit.\n(Remember: < is computed the same as <= and similarly for > and >=)\n")
# getvalidinput is a function that does exactly what it sounds like. See documentation above.
dir <- getvalidinput("Enter your selection: ", "Invalid selection! Choose between 1, 2, 3 or 4!\n", validinput = c('1', '2', '3', '4', 'q'))
if(dir == 'q') { cat("Quitting!\n"); break() }
dir <- as.numeric(dir)
switch(dir,
{
# less than
stat <- getvalidinput("Enter the value for the upper bound c: ", "Test statistic must be a number!\n", validinput='q', low=-1e8, high=1e8);
if(stat == 'q') { cat("Quitting!\n"); break() }
stat <- round(as.numeric(stat), 4);
pv <- round(pnorm(stat), 4); # p-value
max <- max(3,abs(stat)+1); # width for x-axis.
# generate the shaded area.
xcord <- c(-max,seq(-max,stat,.01),stat);
ycord <- c(0,dnorm(seq(-max,stat,.01)),0);
if(pv == 0)
ti <- substitute(paste("P(Z < ", z, ") < .0001", sep=''), list(z=stat)) # create title.
else if(pv == 1)
ti <- substitute(paste("P(Z < ", z, ") > .9999 ", sep=''), list(z=stat)) # create title.
else
ti <- substitute(paste("P(Z < ", z, ") = ", p, sep=''), list(z=stat, p=pv)) # create title.
},
{
# greater than
stat <- getvalidinput("Enter the value for the lower bound c: ", "Test statistic must be a number!\n", validinput='q', low=-1e8, high=1e8);
if(stat == 'q') { cat("Quitting!\n"); break("Quitting!\n") }
stat <- round(as.numeric(stat), 4);
pv <- round(1 - pnorm(stat), 4);
max <- max(3,abs(stat)+1);
xcord <- c(stat,seq(stat,max,.01),max);
ycord <- c(0,dnorm(seq(stat,max,.01)),0);
if(pv == 0)
ti <- substitute(paste("P(Z > ", z, ") < .0001", sep=''), list(z=stat)) # create title.
else if(pv == 1)
ti <- substitute(paste("P(Z > ", z, ") > .9999 ", sep=''), list(z=stat)) # create title.
else
ti <- substitute(paste("P(Z > ", z, ") = ", p, sep=''), list(z=stat, p=pv)) # create title.
},
{
# between
stat1 <- getvalidinput("First, enter the value for the lower bound a:" , "Test statistic must be a number!\n", validinput='q', low=-1e8, high=1e8)
if(stat1 == 'q') { cat("Quitting!\n"); break() }
stat1 <- round(as.numeric(stat1), 4);
cat("Next, enter a value for the upper bound b (must be greater than a).\n");
while(TRUE)
{
stat2 <- getvalidinput("Upper bound: " , "Test statistic must be a number!\n", validinput='q', low=-1e8, high=1e8)
if(stat2 == 'q') { cat("Quitting!\n"); break() }
stat2 <- round(as.numeric(stat2), 4);
if(stat1 < stat2) { break() } else { cat(paste("Upper bound must be greater than lower bound of ", stat1, "!\n", sep='')) }
}
if(stat2 == 'q') { break() } #need second before first only breaks out of inner while
pv <- round(pnorm(stat2) - pnorm(stat1), 4);
max <- max(3,abs(stat1)+1,abs(stat2)+1);
xcord <- c(stat1,seq(stat1,stat2,.01),stat2);
ycord <- c(0,dnorm(seq(stat1,stat2,.01)),0)
if(pv == 0)
ti <- substitute(paste("P(", z1, "< Z < ", z2, ") < .0001", sep=''), list(z1=stat1, z2=stat2))
else if(pv == 1)
ti <- substitute(paste("P(", z1, "< Z < ", z2, ") > .9999 ", sep=''), list(z1=stat1, z2=stat2))
else
ti <- substitute(paste("P(", z1, "< Z < ", z2, ") = ", p, sep=''), list(z1=stat1, z2=stat2, p=pv))
},
{
# outside
stat1 <- getvalidinput("First, enter the value for the lower bound a:" , "Test statistic must be a number!\n", validinput='q', low=-1e8, high=1e8)
if(stat1 == 'q') { cat("Quitting!\n"); break() }
stat1 <- round(as.numeric(stat1), 4);
cat("Next, enter a value for the upper bound b (must be greater than a).\n");
bounds_ok <- FALSE;
while(!bounds_ok)
{
stat2 <- getvalidinput("Upper bound: " , "Test statistic must be a number!\n", validinput='q', low=-1e8, high=1e8)
if(stat2 == 'q') { cat("Quitting!\n"); break() }
stat2 <- round(as.numeric(stat2), 4);
if(stat1 < stat2) { break() } else { cat(paste("Upper bound must be greater than lower bound of ", stat1, "!\n", sep='')) }
}
if(stat2 == 'q') { break() }
pv <- round(1 - pnorm(stat2) + pnorm(stat1), 4);
max <- max(3,abs(stat1)+1,abs(stat2)+1);
xcord1 <- c(-max,seq(-max,stat1,.01),stat1);
ycord1 <- c(0,dnorm(seq(-max,stat1,.01)),0)
xcord2 <- c(stat2,seq(stat2,max,.01),max);
ycord2 <- c(0,dnorm(seq(stat2,max,.01)),0)
ti <- substitute(paste("P(Z < ", z1, " or Z >", z2, ") = ", p, sep=''), list(z1=stat1, z2=stat2, p=pv))
if(pv == 0)
ti <- substitute(paste("P(Z < ", z1, " or Z >", z2, ") < .0001", sep=''), list(z1=stat1, z2=stat2))
else if(pv == 1)
ti <- substitute(paste("P(Z < ", z1, " or Z >", z2, ") > .9999 ", sep=''), list(z1=stat1, z2=stat2))
else
ti <- substitute(paste("P(Z < ", z1, " or Z >", z2, ") = ", p, sep=''), list(z1=stat1, z2=stat2, p=pv))
})
bound <- floor(max*10)/10; # bounds for axis, "floored" to nearest .1.
plot(function(x) dnorm(x),-max,max,
ylim=c(0,.4),
axes=F,
xlab="Z values",
ylab='Density',
main="N(0,1) Distribution",
xaxs='i', #to make the axes start exactly at 0,0
yaxs='i')
axis(1,at=-ceiling(bound):ceiling(bound))
switch(dir,
{
axis(1,at=stat,lwd.ticks=2,col.ticks='red',col.axis='red', font=2,lwd=0, tck=-.04, padj=1.25);
polygon(xcord,ycord,col='red')
},
{
axis(1,at=stat,lwd.ticks=2,col.ticks='red',col.axis='red', font=2,lwd=0, tck=-.04, padj=1.25);
polygon(xcord,ycord,col='red')
},
{
axis(1,at=c(stat1, stat2),lwd.ticks=2,col.ticks='red',col.axis='red', font=2,lwd=0, tck=-.04, padj=1.25);
polygon(xcord,ycord,col='red')
},
{
axis(1,at=c(stat1, stat2),lwd.ticks=2,col.ticks='red',col.axis='red', font=2,lwd=0, tck=-.04, padj=1.25);
polygon(xcord1,ycord1,col='red')
polygon(xcord2,ycord2,col='red')
})
axis(2,at=c(0,1))
legend('topright', eval(ti), col='red', pch=15, bty='n', pt.cex=2)
saveplots <- getvalidinput("Do you want to save your plot to the desktop? y/n: ", "Enter (y)es or (n)o!\n", validinput=c('y','n', 'q'))
if(saveplots == 'q') { cat("Quitting!\n"); break() }
if(saveplots == 'y')
{
titlename <- readline("Please enter your name for the plot: ")
while(TRUE)
{
if(nchar(titlename) > 0) break()
titlename <- readline("Must enter a name! Enter your name: ")
}
os <- .Platform$OS
td <- gsub(":", ".", as.character(Sys.time()))
td <- gsub(" ", "_", td)
if(os == 'unix') #unix == linux or mac
{
name <- paste("~/Desktop/probability_Rplot_", td, sep="")
}
if(os == 'windows')
{
name <- strsplit(as.character(Sys.getenv('HOME')), 'Documents')
name <- paste(name, "Desktop\\probability_Rplot_", td, sep="")
}
name <- paste(name, ".jpg", sep="")
jpeg(name)
plot(function(x) dnorm(x),-max,max,
ylim=c(0,.4),
axes=F,
xlab="Z values",
ylab='Density',
main=paste("N(0,1) Distribution",titlename, sep='\nBy '),
xaxs='i', #to make the axes start exactly at 0,0
yaxs='i')
axis(1,at=-ceiling(bound):ceiling(bound))
switch(dir,
{
axis(1,at=stat,lwd.ticks=2,col.ticks='red',col.axis='red', font=2,lwd=0, tck=-.04, padj=1.25);
polygon(xcord,ycord,col='red')
},
{
axis(1,at=stat,lwd.ticks=2,col.ticks='red',col.axis='red', font=2,lwd=0, tck=-.04, padj=1.25);
polygon(xcord,ycord,col='red')
},
{
axis(1,at=c(stat1, stat2),lwd.ticks=2,col.ticks='red',col.axis='red', font=2,lwd=0, tck=-.04, padj=1.25);
polygon(xcord,ycord,col='red')
},
{
axis(1,at=c(stat1, stat2),lwd.ticks=2,col.ticks='red',col.axis='red', font=2,lwd=0, tck=-.04, padj=1.25);
polygon(xcord1,ycord1,col='red')
polygon(xcord2,ycord2,col='red')
})
axis(2,at=c(0,1))
legend('topright', eval(ti), col='red', pch=15, bty='n', pt.cex=2)
dev.off()
}
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.