# graph a skewed sampling distribution
gskew <- function(n = 1, mu = 0, sd = 1, skew = 0)
{
	se <- sd/sqrt(n)
	x <- seq(mu - 4 * se, mu + 4 * se, length = 501)
	y <- dskew(x, n, mu, sd, skew)
	plot(x, y, xlab = "", ylab = "density", 
		main = paste("Distribution of x-bar, sample size =",n), type = "n")
	lines(x, y, type = "l", col = 4)
	abline(h = 0)
	lines(x, dnorm(x, mu, se), col = 8)
	return(invisible(y))
}

# density of skewed distribution
dskew <- function(x, n = 1, mu = 0, sd = 1, skew = 0)
{
	se <- sd/sqrt(n)
	if(skew == 0)
		y <- dnorm(x, mu, se)
	else {
		a <- (skew * sd)/2
		b <- mu - (2 * sd)/skew
		alpha <- 4/skew/skew
		if(skew > 0)
			y <- (n * dgamma((n * (x - b))/a, n * alpha))/a
		else y <- ( - n * dgamma((n * (x - b))/a, n * alpha))/a
	}
	return(y)
}

# density of symmetric bimodal distribution
dbimod <- function(x,n=1,mu=0,sd=1,d=1)
{
	if(d >= sd)
		stop("sd must be larger than d")
	temp <- matrix(0,nrow=n+1,ncol=length(x))
	for(i in 0:n){temp[i+1,] <- dbinom(i,n,.5)*dnorm(x,mu-d+2*i*d/n,sqrt((sd*sd-d*d)/n))}
	return(apply(temp,2,sum))
}

# graph sampling distribution of symmetric bimodal distribution
gbimod <- function(n=1, mu=0, sd=1, d=0.5)
{
	if(d >= sd)
		stop("sd must be larger than d")
	se <- sd/sqrt(n)
	x <- seq(mu - 4 * se, mu + 4 * se, length = 201)
	y <- dbimod(x, n, mu, sd, d)
	xx <- c(mu - 4 * se, mu + 4 * se)
	yy <- c(0, max(y, dnorm(mu, mu, se)))
	plot(xx, yy, xlab = "", ylab = "density", 
		main = paste("Distribution of x-bar, sample size =",n), type = "n")
	lines(x, y, type = "l", col = 4)
	abline(h = 0)
	lines(x, dnorm(x, mu, se), col = 8)
	return(invisible(y))
}

# graph binomial distribution
gbinom <- function(n, p, low=0, high=n, scale = F)
{
  sd <- sqrt(n * p * (1 - p))
  if(scale && (n > 10)) {
	low <- max(0, round(n * p - 4 * sd))
	high <- min(n, round(n * p + 4 * sd))
  }
  values <- low:high
  probs <- dbinom(values, n, p)
  plot(c(low,high), c(0,max(probs)), type = "n", xlab = "Possible Values", ylab = 
    "Probability", main = paste("Binomial Distribution \n", "n =", 
    n, ", p =", p))
  lines(values, probs, type = "h", col = 2)
  abline(h=0,col=3)
  return(invisible())
}