.packageName <- "Bolstad"
binobp<-function(x,n,a = 1, b = 1 , ret = FALSE){

	# n - the number of trials in the binomial
	# x - the number of observed successes
	# a,b  - the parameters of the Beta prior density (must be > 0)
	# ret - if true then the prior, likelihood, posterior, mean, variance and
  # std. deviation are returned as a list

	if(x>n)
		stop("The number of observed successes (x) must be smaller than the number of trials (n)")
	if(a<=0||b<=0)
		stop("The parameters of the prior must be greater than zero")

	theta<-seq(0.01,0.999,by=0.001)
	prior<-dbeta(theta,a,b)
	likelihood<-dbinom(x,n,prob=theta)
	posterior<-dbeta(theta,a+x,b+n-x)

	plot(theta,posterior,ylim=c(0,1.1*max(posterior,prior)),type="l"
		,lty=1
		,xlab=expression(theta)
		,ylab="Density")
	lines(theta,prior,lty=2)
	
	left<-min(theta)+diff(range(theta))*0.05
	legend(left,max(posterior,prior),lty=1:2,legend=c("Posterior","Prior"))
	
	m1<-(a+x)/(a+b+n)
	v1<-m1*(1-m1)/(a+b+n+1)
	s1<-sqrt(v1)

	cat(paste("Posterior Mean           : ",round(m1,7),"\n"))
	cat(paste("Posterior Variance       : ",round(v1,7),"\n"))
	cat(paste("Posterior Std. Deviation : ",round(s1,7),"\n"))

	probs<-c(0.005,0.01,0.025,0.05,0.5,0.95,0.975,0.99,0.995)
	qtls<-qbeta(probs,a+x,b+n-x)
	names(qtls)<-probs

	cat("\nProb.\tQuantile \n")
	cat("\------\t---------\n")
	for(i in 1:length(probs))
		cat(paste(round(probs[i],3),"\t",round(qtls[i],7),"\n",sep=""))


	if(ret)
		return(list(posterior=posterior,likelihood=likelihood,prior=prior,theta=theta,mean=m1,var=v1,sd=s1,quantiles=qtls))
	
}
binodp<-function(x,n,uniform = TRUE, n.theta = 10, theta = NULL, theta.prior = NULL, ret = FALSE){

	# n - the number of trials in the binomial
	# x - the number of observed successes
	# theta - the probability of success 
	# theta.prior - the associated prior probability mass
	# ret - if true then the likelihood and posterior are returned as a 
	# list
	
	if(x>n)
		stop("The number of observed successes (x) must be smaller than the number of trials (n)") 
	if(n.theta<3)
		stop("Number of prior values of theta must be greater than 2")

	if(is.null(theta)&is.null(theta.prior)&uniform){
		theta<-seq(0,1, length = n.theta)
		theta.prior<-rep(1/n.theta,n.theta)
	}

	if(sum(theta<0) > 0 | sum(theta > 1) > 0) # check that probabilities lie on [0,1]
		stop("Values of theta must be between 0 and 1 inclusive")

	if(sum(theta.prior<0)>0|sum(theta.prior>1)>1)
		stop("Prior probabilities must be between 0 and 1 inclusive")

	if(round(sum(theta.prior),7)!=1){
		warning("The prior probabilities did not sum to 1, therefore the prior has been normalized")
		theta.prior<-theta.prior/sum(theta.prior)
	}

	if(!uniform&(is.null(theta)|is.null(theta.prior)))
		stop("If you wish to use a non-uniform discrete prior then you must supply a theta vector and an associated probability vector theta.prior")


	n.theta<-length(theta)
	likelihood<-dbinom(x,n,theta)*theta.prior


	posterior<-likelihood/sum(likelihood)

	plot(theta,posterior,ylim=c(0,1.1*max(posterior,theta.prior)),pch="o",
		xlab=expression(theta),ylab=expression(Probabilty(theta)))
	points(theta,theta.prior,pch="+")

	legend(max(c(0.05,min(theta))),max(posterior,theta.prior),pch=c("o","+"),legend=c("Posterior","Prior"))

	# calculate the Conditional distribution

	f.cond<-matrix(0,nrow=n.theta,ncol=n+1)
	rownames(f.cond)<-as.character(round(theta,3))
	colnames(f.cond)<-as.character(0:n)

	for(i in 1:n.theta)
		f.cond[i,]<-dbinom(0:n,n,theta[i])

	cat("Conditional distribution of x given theta and  n:\n\n")
	print(round(f.cond,4))	

	# caculate the joint distribution of theta and x given n

	f.joint<-diag(theta.prior)%*%f.cond
	cat("\nJoint distribution:\n\n")
	print(round(f.joint,4))	

	# calculate the marginal distribtion

	f.marg<-matrix(1,nrow=1,ncol=n.theta)%*%f.joint
	cat("\nMarginal distribution of x:\n\n")
	print(round(f.marg,4))
	cat("\n\n")	

	# finally display the prior, likelihood, and posterior

	results<-cbind(theta.prior,likelihood,posterior)
	rownames(results)<-as.character(round(theta,3))
	colnames(results)<-c("Prior","Likelihood","Posterior")

	print(results)






	if(ret)
		return(list(f.cond=f.cond,f.joint=f.joint,f.marg=f.marg,likelihood=likelihood,posterior=posterior,theta=theta,theta.prior=theta.prior))
	
}
binogcp<-function(x,n,density="uniform", params = c(0,1), n.theta = 1000, theta = NULL, theta.prior = NULL, ret = FALSE){

	# n - the number of trials in the binomial
	# x - the number of observed successes
	# density - may be one of "exp","normal","uniform" or "user"
	# params - if the density is not "user" then a vector of parameters
	# must be supplied. 
	#	exp:		rate
	#	normal: 	mean,sd
	#	uniform: 	min,max
	# n.theta - the number of points to divide the [0,1] interval into  

	# theta and theta.prior are only specified if density == "user"
	# theta - the probability of success 
	# theta.prior - the associated prior probability mass
	# ret - if true then the likelihood and posterior are returned as a 
	# list
	
	if(x>n)
		stop("The number of observed successes (x) must be smaller than the number of trials") 
	if(n.theta<100)
		stop("Number of prior values of theta must be greater than 100")
	
	if(is.null(theta)||is.null(theta.prior))
		theta<-seq(0+1/n.theta,1-1/n.theta,length=n.theta)
	else{
		if(length(theta)!=length(theta.prior))
			stop("theta and theta.prior must have same length")
		
		if(sum(theta<0|theta>1)>0) # check that probabilities lie on [0,1]
			stop("Values of theta must be between 0 and 1 inclusive")
	}

	if(density=="beta"){
		if(length(params)<2){
			warning("Beta prior requires two shape parameters. Default value Beta(1,1) = Uniform is being used")
			a<-1
			b<-1
		}else{
			if(params[1]<=0|params[2]<0)
				stop("Beta prior shape parameters must be greater than zero")
			a<-params[1]
			b<-params[2]
		}
		theta.prior<-dbeta(theta,a,b)
	}else	if(density=="exp"){
		if(params[1]<=0){
			stop("Parameter for exponential density must be greater than zero")
		}else{
			rate<-params[1]
			theta.prior<-dexp(theta,rate)
		}
	}else if(density=="normal"){
		if(length(params)<2)
			stop("Normal prior requires a mean and std. deviation")
		else{
			mx<-params[1]
			sx<-params[2]
			if(sx<=0)
				stop("Std. deviation for normal prior must be greater than zero")
			theta.prior<-dnorm(theta,mx,sx)
		}
	}else if(density=="uniform"){
		if(length(params)<2)
			stop("Uniform prior requires a minimum and a maximum")
		else{
			minx<-params[1]
			maxx<-params[2]
			
			if(maxx<=minx)
				stop("Maximum must be greater than minimum for a uniform prior")
			theta.prior<-dunif(theta,minx,maxx)
		}
	}else if (density!="user"){
		stop(paste("Unrecognized density :",density))
	}

	likelihood<-(theta^x)*((1-theta)^(n-x))

	# Numerically integrate the denominator
	# First calculate the height of the function to be integrated

	f.x.theta<-likelihood*theta.prior

	# Now get a linear approximation so that we don't have to worry about
	# the number of points specified by the user

	ap<-approx(theta,f.x.theta,n=513)
	integral<-sum(ap$y[2*(1:256)-1]+4*ap$y[2*(1:256)]+ap$y[2*(1:256)+1])
	integral<-(ap$x[2]-ap$x[1])*integral/3

	posterior<-likelihood*theta.prior/integral

	plot(theta,posterior,ylim=c(0,1.1*max(posterior,theta.prior)),lty=1,type="l",col="blue",
		xlab=expression(theta),ylab="Density")
	lines(theta,theta.prior,lty=2,col="red")
	
	left<-min(theta)+diff(range(theta))*0.05
	legend(left,max(posterior,theta.prior),lty=1:2,col=c("blue","red"),legend=c("Posterior","Prior"))
	if(ret)
		return(list(likelihood=likelihood,posterior=posterior,theta=theta,theta.prior=theta.prior))
	
}
normdp<-function(x,sigma.x,uniform = TRUE, n.mu = 50, mu = NULL, mu.prior = NULL, ret = FALSE){

	# x - the vector of observations
	# sigma.x - the population standard deviation
	# uniform - default to a discrete uniform prior
	# mu - vector of possible values of the population mean 
	# mu.prior - the associated prior probability mass
	# ret - if true then the likelihood and posterior are returned as a 
	# list
	
	if(n.mu<3)
		stop("Number of prior values of theta must be greater than 2")

	if(uniform){
		mu<-seq(min(x)-sigma.x,max(x)+sigma.x,length = n.mu)
		mu.prior<-rep(1/n.mu,n.mu)
	}

	if(any(mu.prior<0) | any(mu.prior>1))
		stop("Prior probabilities must be between 0 and 1 inclusive")

	if(round(sum(mu.prior),7)!=1){
		warning("The prior probabilities did not sum to 1, therefore the prior has been normalized")
		mu.prior<-mu.prior/sum(mu.prior)
	}

	if(!uniform&(is.null(mu)|is.null(mu.prior)))
		stop("If you wish to use a non-uniform discrete prior then you must supply a mean vector, mu, and an associated probability vector, mu.prior")


	n.mu<-length(mu)
	mx<-mean(x)
	nx<-length(x)
	snx<-sigma.x^2/nx
	likelihood<-exp(-0.5*(mx-mu)^2/snx)

	posterior<-likelihood*mu.prior/sum(likelihood*mu.prior)

	plot(mu,posterior,ylim=c(0,1.1*max(posterior,mu.prior)),pch="o",
		xlab=expression(mu),ylab=expression(Probabilty(mu)))
	points(mu,mu.prior,pch="+")
	
	left<-min(mu)+diff(range(mu))*0.05
	legend(left,max(posterior,mu.prior),pch=c("o","+"),legend=c("Posterior","Prior"))

	if(ret)
		return(list(likelihood=likelihood,posterior=posterior,mu=mu,mu.prior=mu.prior))
	
}
normgcp<-function(x,sigma.x,density="uniform" , params = NULL, n.mu = 50, mu = NULL, mu.prior = NULL, ret = FALSE){

	# x - the vector of observations
	# sigma.x - the population standard deviation
	# density - distributional form of the prior density
	# can be one of : normal, unform, or user 
	# by default a continuous uniform prior is used
	# mu - vector of possible values of the population mean 
	# mu.prior - the associated prior probability mass
	# ret - if true then the likelihood and posterior are returned as a 
	# list

	mean.x<-mean(x)
	
	if(n.mu<3)
		stop("Number of prior values of mu must be greater than 2")

	if(density=="normal"){
		if(is.null(params)|length(params)<1)
			stop("You must supply a mean for a normal prior")
		mx<-params[1]

		if(length(params)==2) # user has supplied sd as well
			s.x<-params[2]
		else
			s.x<-sigma.x

		mu<-seq(mx-3.5*s.x,mx+3.5*s.x,length = n.mu)
		mu.prior<-dnorm(mu,mx,s.x)
	} else if(density=="uniform"){
		if(is.null(params)){
			# set params to mean+/-3.5sd by default
			params<-c(mean.x-3.5*sigma.x,mean.x+3.5*sigma.x)
		}
		if(length(params)<2)
			stop("You must supply a minimum and a maximum to use a uniform prior")
		minx<-params[1]
		maxx<-params[2]
		if(maxx<=minx)
			stop("The maximum must be greater than the minimum for a uniform prior")
		mu<-seq(minx,maxx,length = n.mu)
		mu.prior<-dunif(mu,minx,maxx)
	}else{
		# user specified prior
		if(is.null(mu)|is.null(mu.prior))
			stop("If you wish to use a non-uniform continuous prior then you must supply a mean vector, mu, and an associated density vector, mu.prior")
	}

	if(any(mu.prior<0) | any(mu.prior>1))
		stop("Prior probabilities must be between 0 and 1 inclusive")

	crude.int<-sum(diff(mu)*mu.prior[-1])
	if(round(crude.int,3)!=1){
		warning("The prior probabilities did not sum to 1, therefore the prior has been normalized")
		mu.prior<-mu.prior/crude.int
		print(crude.int)
	}

	n.mu<-length(mu)
	mx<-mean(x)
	nx<-length(x)
	snx<-sigma.x^2/nx
	likelihood<-exp(-0.5*(mx-mu)^2/snx)

	# Numerically integrate the denominator
	# First calculate the height of the function to be integrated

	f.x.mu<-likelihood*mu.prior

	# Now get a linear approximation so that we don't have to worry about
	# the number of points specified by the user

	ap<-approx(mu,f.x.mu,n=513)
	integral<-sum(ap$y[2*(1:256)-1]+4*ap$y[2*(1:256)]+ap$y[2*(1:256)+1])
	integral<-(ap$x[2]-ap$x[1])*integral/3

	posterior<-likelihood*mu.prior/integral

	plot(mu,posterior,ylim=c(0,1.1*max(posterior,mu.prior)),type="l",
		lty=1,col="blue",
		xlab=expression(mu),ylab=expression(Probabilty(mu)))
	lines(mu,mu.prior,lty=2,col="red")

	left<-min(mu)+diff(range(mu))*0.05
	legend(left,max(posterior,mu.prior),lty=1:2,col=c("blue","red"),legend=c("Posterior","Prior"))

	if(ret)
		return(list(likelihood=likelihood,posterior=posterior,mu=mu,mu.prior=mu.prior))
	
}
normnp<-function(x,sigma.x,m.x=0,s.x=1,n.mu = 100,ret=FALSE){

	# x - the vector of observations
	# sigma.x - the population standard deviation
	# m.x - the mean of the normal prior
	# s.x - the standard deviation of the normal prior
	# ret - if true then the prior, likelihood, posterior, mean, variance, and
	# quantiles are returned as a list

	mean.x<-mean(x)
	
	if(n.mu<100)
	{
		warning("Number of prior values of mu must be greater than 100")		
		n.mu<-100
	}

	if(s.x<=0)
		stop("The std. deviation of the prior must be greater than zero")

	lb<-m.x-3.5*s.x
	ub<-m.x+3.5*s.x

	mu<-seq(lb,ub,length=n.mu)
	mu.prior<-dnorm(mu,m.x,s.x)
	
	n.x<-length(x)
	
	likelihood<-exp(-n.x/(2*sigma.x^2)*(mean.x-mu)^2)

	precision<-1/s.x^2
	post.precision<-precision+(n.x/sigma.x^2)
	post.sd<-sqrt(1/post.precision)
	post.mean<-(precision/post.precision*m.x)+((n.x/sigma.x^2)/post.precision*mean.x)

	cat(paste("Posterior mean           : ",round(post.mean,7),"\n",sep=""))
	cat(paste("Posterior std. deviation : ",round(post.sd,7),"\n",sep=""))

	posterior<-dnorm(mu,post.mean,post.sd)
	
	plot(mu,posterior,ylim=c(0,1.1*max(posterior,mu.prior)),type="l",
		lty=1,col="blue",
		xlab=expression(mu),ylab=expression(Probabilty(mu)))
	lines(mu,mu.prior,lty=2,col="red")

	left<-min(mu)+diff(range(mu))*0.05
	legend(left,max(posterior,mu.prior),lty=1:2,col=c("blue","red"),legend=c("Posterior","Prior"))

	probs<-c(0.005,0.01,0.025,0.05,0.5,0.95,0.975,0.99,0.995)
	qtls<-qnorm(probs,post.mean,post.sd)
	names(qtls)<-probs

	cat("\nProb.\tQuantile \n")
	cat("\------\t---------\n")
	for(i in 1:length(probs))
		cat(paste(round(probs[i],3),"\t",round(qtls[i],7),"\n",sep=""))

	if(ret)
		return(list(prior=prior,likelihood=likelihood,posterior=posterior,mean=post.mean,sd=post.sd,quantiles=qtls))
	
}
sintegral<-function(x,fx,n.pts=256,ret=FALSE)
{
	# numerically integrates fx over x using Simpsons rule
	# x - a sequence of x values
	# fx - the value of the function to be integrated at x
	# n.pts - the number of points to be used in the integration
	# ret - if true returns the partial sums of the integration


	n.x<-length(x)

	if(n.x!=length(fx))
		stop("Unequal input vector lengths")

	if(n.pts<64)
		n.pts<-64

	# use linear approximation to get equally spaced x values


	ap<-approx(x,fx,n=2*n.pts+1)

	h<-diff(ap$x)[1]

	integral<-h*(ap$y[2*(1:n.pts)-1]
			+4*ap$y[2*(1:n.pts)]
			+ap$y[2*(1:n.pts)+1])/3

	if(ret)
		return(list(x=ap$x[2*(1:n.pts)],y=cumsum(integral)))

	return(sum(integral))
}
sscsample<-function(size, n.samples, sample.type="simple", x = NULL, 
					strata = NULL, cluster = NULL, cl.size = NULL, ret=FALSE)
{
	# size - the sample size 
	# n.samples - the number of samples to draw
	# sample.type - the method of sampling can be one of :
	#      "cluster","stratified","simple"

	# DO NOT set the following parameters unless you know what you're doing!
	# x - a data vector
	# strata - the stratum each data point in x belongs to
	# cluster - the cluster each data point in x belongs to
	# cl.size - the number of clusters to sample 
	#	- must be less than the number of clusters

	# if ret is true then the samples and their summary statistics are
	# return in a list
	
	data(sscsample.data)

	if(is.null(x))
		x<-sscsample.data$value

	nx<-length(x)

	if(size>nx)
		stop("Sample size must be less than population size")

	if(is.null(strata))
		strata<-sscsample.data$stratum
	
	# just in case the strata numbers are not evenly spaced
	strata.names<-unique(strata)
	n.strata<-length(strata.names)

	if(nx!=length(strata))
		stop("The length of the strata and data vectors must be equal")

	if(is.null(cluster))
		cluster<-sscsample.data$cluster

	n.clusters<-length(unique(cluster))

	if(nx!=length(cluster))
		stop("The length of the cluster and data vectors must be equal")

	samples<-matrix(0,nrow=size,ncol=n.samples)

	for(r in 1:n.samples){
		if(sample.type=="simple"|sample.type==1){
			sample.idx<-sample(1:nx,size)

		} else if(sample.type=="stratified"|sample.type==2){
			
			stratified.data<-split(1:nx,strata.names)

			k<-size*sapply(stratified.data,length)/nx
			
			for(stratum in 1:n.strata){
				if(stratum==1)
					sample.idx<-sample(stratified.data[[stratum]],k[stratum])
				else
					sample.idx<-c(sample.idx,sample(stratified.data[[stratum]],k[stratum]))
			}	
		} else if(sample.type=="cluster"|sample.type==3){
			# just in case the cluster numbers are not evenly spaced			
			
			cluster.names<-unique(cluster) 

			if(cl.size<0)
				cl.size<-4

			if(cl.size>n.clusters)
				stop("The number of randomly sampled clusters must be less that the overall number of clusters")
		
			cluster.sample<-sample(cluster.names,cl.size)

			clustered.data<-split(1:nx,cluster.names)
			
			for(i in 1:cl.size)
			{
				cl<-cluster.sample[i]
				if(i==1)
					sample.idx<-clustered.data[[cl]]
				else
					sample.idx<-c(sample.idx,clustered.data[[cl]])
			}		
		} else
			stop(paste("Unknown sampling sample.type :",sample.type))
	
		samples[,r]<-sample.idx
	}

	means<-rep(0,n.samples)
	s.strata<-matrix(0,nrow=n.samples,ncol=n.strata)
	sample.out<-matrix(0,nrow=size,ncol=n.samples)

	cat("Sample\tMean   \tStratum 1\tStratum 2\tStratum 3\n")
	cat("------\t-------\t---------\t---------\t---------\n")

	for(r in 1:n.samples)
	{
		idx<-samples[,r]
		means[r]<-mean(x[idx])

		for(j in 1:n.strata)
			s.strata[r,j]<-sum(strata[idx]==strata.names[j])

		sample.out[,r]<-x[idx]

		cat(paste(r,"\t",round(means[r],4),"\t",s.strata[r,1],"\t\t",s.strata[r,2]
			,"\t\t",s.strata[r,3],"\n",sep=""))
	}
	
	if(ret)
		return(list(samples=samples,s.strata=s.strata,means=means))
	
}
xdesign<-function(x=NULL,y=NULL,corr=0.8,size=20,n.treatments=4,n.rep=500)
{
	if(is.null(x)) # simulate some data
	{
		nx<-size*n.treatments
		x<-rnorm(nx)
		y<-rnorm(nx)

		y<-sqrt(1-corr^2)*y+corr*x
	}
	nx<-size*n.treatments
	if(length(x)!=length(y))
		stop("x and y must be of equal length")
	if(length(x)!=size*n.treatments)
		stop("x and y must be equal to the same size times the number of treatments")
	if(corr<(-1)|corr>1)
		stop("Correlation coeficient must be between -1 and 1")

	if(n.rep<10)
		stop("Must have at least 10 Monte Carlo replicates")

	cat("Variable\tN\tMean\tMedian\tTrMean\tStDev\tSE Mean\n")
	cat(paste("X\t",length(x),
									round(mean(x),3),
									round(median(x),3),
									round(mean(x,trim=0.1),3),
									round(sd(x),3),
									round(sd(x)/sqrt(length(x)),3),sep="\t"))
	cat("\n")
	cat(paste("Y\t",length(y),
									round(mean(y),3),
									round(median(y),3),
									round(mean(y,trim=0.1),3),
									round(sd(y),3),
									round(sd(y)/sqrt(length(y)),3),sep="\t"))
	cat("\n\n")

	qx<-quantile(x,c(0.25,0.75))
	qy<-quantile(y,c(0.25,0.75))

	cat("Variable\tMinimum\tMaximum\tQ1\tQ3\n")
	cat(paste("X\t",round(min(x),3)
								,round(max(x),3)
								,round(qx[1],3)
								,round(qx[2],3),sep="\t"))
	cat("\n")
	cat(paste("Y\t",round(min(y),3)
								,round(max(y),3)
								,round(qy[1],3)
								,round(qy[2],3),sep="\t"))

	cat("\n\n")

	cat("The Pearson correlation between Y and Y is: ");
	cat(paste(round(cor(x,y),3),"\n\n"))

	plot(x,y)

	ssx<-rep(0,n.rep)
	ssy<-rep(0,n.rep)

	treat.groupmean<-matrix(0,ncol=n.treatments,nrow=n.rep)
	block.groupmean<-matrix(0,ncol=n.treatments,nrow=n.rep)

	for(block in c(FALSE,TRUE))
	{
		# block is indicator for blocking
		# FALSE =  completely randomized design, 
		# TRUE = randomized block design

		for(i in 1:n.rep)
		{
			if(!block)
			{
				group<-rep(1:n.treatments,size)

				z<-rnorm(nx)
				o<-order(z)
				z<-z[o]
				group<-group[o]

				x2<-x
				y2<-y
			}
			else
			{
				o<-order(x)
				x2<-x[o]
				y2<-y[o]
				
				group<-NULL

				for(j in 1:size)
				{
					gp<-1:n.treatments
					z<-rnorm(n.treatments)
					gp<-gp[order(z)]
					group<-c(group,gp)
				}
			}

			split.x<-split(x2,group)
			split.y<-split(y2,group)

			x.bar<-sapply(split.x,mean)
			y.bar<-sapply(split.y,mean)

			x.mean<-mean(x.bar)
			y.mean<-mean(y.bar)

			ssx[i]<-sum((x.bar-x.mean)^2)
			ssy[i]<-sum((y.bar-y.mean)^2)

			treat.groupmean[i,]<-y.bar
			block.groupmean[i,]<-x.bar
		}
		
		if(!block)
		{
			treat.var0<-as.vector(treat.groupmean)
			block.var0<-as.vector(block.groupmean)
			index0<-rep(1:n.treatments,rep(n.rep,n.treatments))
		}
		else
		{
			treat.var1<-as.vector(treat.groupmean)
			block.var1<-as.vector(block.groupmean)
			index1<-rep(1:n.treatments,rep(n.rep,n.treatments))
		}
	}

	treat.var<-c(treat.var0,treat.var1)
	block.var<-c(block.var0,block.var1)
	index<-c(index0,index1)

	ind<-rep(1:2,c(length(treat.var0),length(treat.var1)))
	ind<-n.treatments*(ind-1)+index

	par(ask=interactive())
	
	rng<-range(block.var)
	y.lims<-max(abs(c(rng[1]-0.1*diff(rng),rng[2]+0.1*diff(rng))))
	y.lims<-c(-y.lims,y.lims)

	boxplot(block.var~ind
		,main="Boxplots of Lurking/Blocking variable group means"
		,sub="Lurking variable in completely randomized design\nBlocking variable in randomized block design",col=rep(c("blue","red"),rep(n.treatments,2))
		,ylim=y.lims)
	legend(n.treatments+0.5,rng[2],legend=c("Completely randomized design","Randomized block design"),fill=c("blue","red"))
	
	rng<-range(treat.var)
	y.lims<-max(abs(c(rng[1]-0.1*diff(rng),rng[2]+0.1*diff(rng))))
	y.lims<-c(-y.lims,y.lims)
	boxplot(treat.var~ind
		,main="Boxplots of treatment group means"
		,col=rep(c("blue","red"),rep(n.treatments,2))
		,ylim=y.lims)
	legend(n.treatments+0.5,rng[2],legend=c("Completely randomized design","Randomized block design"),fill=c("blue","red"))

	x<-treat.var[ind<=n.treatments]
	y<-treat.var[ind>n.treatments]
	cat("Variable\tN\tMean\tMedian\tTrMean\tStDev\tSE Mean\n")
	cat(paste("Randomized",length(x),
									round(mean(x),3),
									round(median(x),3),
									round(mean(x,trim=0.1),3),
									round(sd(x),3),
									round(sd(x)/sqrt(length(x)),3),sep="\t"))
	cat("\n")
	cat(paste("Blocked\t",length(y),
									round(mean(y),3),
									round(median(y),3),
									round(mean(y,trim=0.1),3),
									round(sd(y),3),
									round(sd(y)/sqrt(length(y)),3),sep="\t"))
	cat("\n\n")

	qx<-quantile(x,c(0.25,0.75))
	qy<-quantile(y,c(0.25,0.75))

	cat("Variable\tMinimum\tMaximum\tQ1\tQ3\n")
	cat(paste("Randomized",round(min(x),3)
								,round(max(x),3)
								,round(qx[1],3)
								,round(qx[2],3),sep="\t"))
	cat("\n")
	cat(paste("Blocked\t",round(min(y),3)
								,round(max(y),3)
								,round(qy[1],3)
								,round(qy[2],3),sep="\t"))

	cat("\n\n")

	return(invisible(list(block.means=block.var,treat.means=treat.var,ind=ind)))
}
