Nothing
#' Illustration of the Normal approximation to Binomial
#'
#' When certain conditions are met (see Details), the Binomial distribution can be approximated by the Normal one. The function \code{AproxBinomNorm} illustrates this fact by plotting the mass diagram corresponding with the discrete distribution (parameters are given by the user) on which the associated Normal density function is also displayed.
#'
#' @param n number of independent Bernouilli trials.
#' @param p probability of success associated with the Bernouilli trial.
#' @param legend logical argument indicating whether to display the legend on the plot or not, default to TRUE.
#' @param xlab x-axis label; default to empty.
#' @param ylab y-axis label; default to "Probability."
#' @param main title; default to "Normal approximation to Binomial".
#' @param col.fill colour to fill-in the bars; default to grey.
#' @param col.line colour to draw the line of the Normal density; default to red.
#' @param lwd line width for the Normal density, a positive number; default to 2.
#'
#' @details The approximation is accurate only if one of these three conditions is met:
#'
#' - p in (0.1, 0.9) and n>=30,
#'
#' - p in [0,0.1] and np>5,
#'
#' - p in [0.9,1] and n(1-p)>5.
#'
#' @return This function is called for the side effect of drawing the plot.
#'
#' @examples
#' n=45; p=0.4
#' AproxBinomNorm(n,p)
#' AproxBinomNorm(n,p,col.fill="blue",col.line="orange")
#' AproxBinomNorm(n,p,legend=FALSE)
#'
#' @export
AproxBinomNorm<-function(n,p,legend=TRUE,xlab="", ylab="Probability", main="Normal approximation to Binomial", col.fill="grey", col.line="red",lwd=2){
if(length(n)!=1){stop("'n' must be a single positive integer")}
if (any(!is.finite(n))){stop("'n' must be a single positive integer")}
if(!is.wholenumber(n)){stop("'n' must be a positive integer")}
if(n<=0){stop("'n' must be a single positive integer")}
if(length(p)!=1){stop("'p' must be a single number in (0,1)")}
if (any(!is.finite(n))){stop("'p' must be a single number in (0,1)")}
if(p<=0 | p>=1){stop("'p' must be a single number in (0,1)")}
if(is.null(legend)){legend=TRUE}
if(!is.null(legend)){if(!is.logical(legend)){stop("'legend' must be a single logical value")}}
if(length(col.line)!=1) stop("The argument 'col.line' must be a single colour")
if(col.line%in%c(NA,NaN,Inf,-Inf)) stop("The argument 'col.line' must be a single colour")
if(length(col.fill)!=1) stop("The argument 'col.fill' must be a single colour")
if(col.fill%in%c(NA,NaN,Inf,-Inf)) stop("The argument 'col.fill' must be a single colour")
if(p>0.1 & p<0.9){if(n<30){warning("The aproximation might not be accurate because p is in (0.1,0.9) but n < 30")}}
if(p<=0.1){if(n*p<=5){warning("The aproximation might not be accurate because p < 0.1 but np <= 5")}}
if(p>=0.9){if(n*(1-p)<=5){warning("The aproximation might not be accurate because p > 0.9 but n(1-p) <= 5")}}
if(is.null(xlab)){xlab=""}
if(is.null(ylab)){ylab="Probability"}
if(is.null(main)){main="Normal approximation to Binomial"}
if(is.null(col.fill)){col.fill="grey"}
if(is.null(col.line)){col.line="red"}
if(is.null(lwd)){lwd=2}
if(!is.numeric(lwd)|length(lwd)!=1) stop("The argument 'lwd' must be a positive integer")
if(!is.finite(lwd)|lwd<=0|lwd!=round(lwd,0)) stop("The argument 'lwd' must be a positive integer")
mu=n*p
sigma=sqrt(n*p*(1-p))
xi=0:n
pi=dbinom(xi,n,p)
#barplot(as.table(pi),names.arg=xi,ylim=c(0,max(pi)+0.03),xlab="",ylab="Probability",main="Normal approximation to binomial")
plot(0:(n+1),type="n", xlim=c(-1,n+1),ylim=c(0,max(pi)+0.05),xlab=xlab,ylab=ylab,axes=F,main=main)
axis(1,at=c(-1,xi,max(xi)+1))
axis(2,at=round(seq(0,(min(1,max(pi)+0.02)),len=5),2))
segments(x0=xi-0.3,y0=0,x1=xi-0.3,y1=pi)
segments(x0=xi+0.3,y0=0,x1=xi+0.3,y1=pi)
segments(x0=xi-0.3, y0=pi, x1=xi+0.3, y1=pi)
points(xi,pi,pch=16)
for(i in 1:(n+1)){polygon(x=c(i-1-0.3,seq(i-1-0.3,i-1+0.3,len=100),i-1+0.3), y=c(0,rep(pi[i],100),0), col=col.fill)}
lines(seq(-1,n,len=1000),dnorm(seq(-1,n,len=1000),mean=mu,sd=sigma),col=col.line,lwd=lwd)
if(legend==TRUE){legend("topright",legend = c(paste("Bi ( ",n," , ",round(p,2)," )",sep=""),paste("N ( ",round(mu,2)," , ",round(sigma,2)," )",sep="")), col=c(col.fill,col.line), lwd=c(NA,lwd), pch=c(15,NA), bty="n")}
# cat(paste("The distributions plotted are Bi(",n,",",round(p,2),") and N(",round(mu,2),",",round(sigma,2),").", sep=""))
}
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.