test_standardDev.R

################################################################################################
# AUTHOR: JD Long
# PURPOSE: put confidence bands around a sample standard deviation
################################################################################################


#based on example code from here: 
#http://www.boost.org/doc/libs/1_36_0/libs/math/doc/sf_and_dist/html/math_toolkit/dist/stat_tut/weg/cs_eg/chi_sq_intervals.html

#thanks to Vince Buffallo for sending me the link and pointing me in the right direction. 
#the sdevConfIntervals takes the following inputs:
#    sample Variance, sampleN, and confidence level (e.g. .95 or .90) 
#    returns the upper and lower bounds for the standard deviation


sdevConfIntervals <- function(sampleVariance, sampleN, confidence){
	upperBound <- sqrt( ((sampleN-1) * sampleVariance) / 
					qchisq((1-confidence)/2, sampleN-1)
	)
	lowerBound <- sqrt( ((sampleN-1) * sampleVariance) / 
					qchisq(1-(1-confidence)/2, sampleN-1)
	)
	return(list(upper=upperBound, lower=lowerBound))	
}



################################################################################################
# the next bit is just some testing code so I could prove to myself it works
################################################################################################

testStdevDist <- function(sampleSize, confidence){
	sample <- rnorm(sampleSize)
	sampleVariance <- var(sample)
	sampleN <- length(sample)
	confidence <- .95
	
	return(sdevConfIntervals(sampleVariance, sampleN, confidence))
}



#test this on 10,000 samples of 100 each 
outData <- NULL
for (i in 1:1e4){
	outData[[i]]<- testStdevDist(100, .9)
}

#test if the actual sdev is between the upper and lower bounds. 
testOut <- NULL
for (i in 1:length(outData)){
	testOut[i] <- if(outData[[i]]$lower <= 1 & outData[[i]]$upper >= 1) {1} else {0}
}

#if everything goes as planned this should be ~ the confidence number (.90, .95, etc)
mean(testOut)