# R/prob.R In learnstats: An Interactive Environment for Learning Statistics

#### Documented in prob

```#' Run the Probability Walkthrough script.
#'
#' @export
#' @examples
#' \dontrun{prob()}
prob<-function()
{
### Probability script for R.
### by Josh Errickson
### 1.4   - Changed title to always N(0,1), give prob in legend.
### 1.3.4 - Never display probabilies as 0 or 1... just say <.0001 or >.9999
### 1.3.3 - Display instructions when script is loaded.
### 1.3.2 - Added name to plots.
### 1.3.1 - Fixes to user input handling.
### 1.3   - Helper function to handle user input - should no longer ever crash! Just keep requesting new input.
### 1.2.1 - Added better error checking for first choice.
### 1.2   - Auto-saving
### 1.1   - Added looping action
### 1.0.2 - Minor typo (double lower bound)
### 1.0.1 - Clarified some language in prompts.
### 1.0   - First version as cribbed from pval.R

##' Given a prompt and acceptable inputs, query user for input until they enter an acceptable value.
##'
##'
##' @title
##' @param prompt String shown to users to prompt for input.
##' @param error Error string shown when user enters an invalid selection.
##' @param validinput If non-NULL, only accept input in this list.
##' @param low Only accept numeric values above this value.
##' @param high Only accept numeric values below this value. If low>=high,
##' do not accept any numeric input not explictly in validinput.
##' @return Acceptable user input.
##' @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)
{
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.

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,
{
polygon(xcord,ycord,col='red')
},
{
polygon(xcord,ycord,col='red')
},
{
polygon(xcord,ycord,col='red')
},
{
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')
{
while(TRUE)
{
if(nchar(titlename) > 0) break()
}
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,
{
polygon(xcord,ycord,col='red')
},
{
polygon(xcord,ycord,col='red')
},
{
polygon(xcord,ycord,col='red')
},
{
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()
}

}
}
```

## Try the learnstats package in your browser

Any scripts or data that you put into this service are public.

learnstats documentation built on May 29, 2017, 9:08 a.m.