probability or statistics - How do I calculate Standard Deviation confidence intervals in Mathematica?
I have a series of simulation results and I would like to demonstrate how the confidence interval for the mean and standard deviation decrease as the number of results increases.
I have done this for the mean using the inbuilt MeanCI function as shown below (I've replaced my simulation results with random data).
Needs["HypothesisTesting`"]
data = Table[RandomReal[{-1, 1}], {i, 1, 100}];
movingData = Table[data[[1 ;; i]], {i, 2, Length[data]}];
meanCI = MeanCI[#, ConfidenceLevel -> .80] & /@ movingData;
ListPlot[Transpose[meanCI]]

What I would like to do is create a similar chart showing how the standard deviation confidence interval decreases. There isn't a standard function for StandardDeviationCI but is there a way to use the tools inside the HypothesisTesting package to get a similar result?
Answer
The standard deviation confidence interval may be obtained via Sqrt@VarianceCI[sample]. If $X_1, X_2,\dots,X_n$ are iid observations from a normal distribution, then $${(n-1)\,s^2\over\sigma^2} \approx \chi^2_{n-1}\,.$$ The confidence interval for $\sigma$ or $\sigma^2$ can be computed from the sample and the $\chi^2$ distribution. This is what VarianceCI does. The sample means of a normal distribution are distributed according to a Student $t$ distribution if $\sigma$ is unknown. The confidence interval can be computed from the $t$ distribution, and this is what MeanCI does.
The OP's example data is from a uniform distribution, not a normal distribution. When the sample size is over 50, the error appears to be practically negligible (compared to the appropriate UniformSumDistribution). I don't know how to compute the exact distribution of the standard deviation, but since the sampling distribution approaches a normal distribution, I will assume that VarianceCI will be adequately approximate for larger sample sizes.
The plot below shows the sample standard deviation $s$ bracketed by the confidence interval for $\sigma$. One may notice that the confidence interval is not symmetric. That is because the \chi^2 distribution is skewed.
Needs["HypothesisTesting`"];
SeedRandom[1];
data = RandomReal[{-1, 1}, 100];
sdci = Table[
Sqrt@VarianceCI[data[[1 ;; i]], ConfidenceLevel -> .80], {i, 2, Length[data]}];
ListPlot[Transpose[sdci] ~Append~
Table[StandardDeviation[data[[1 ;; i]]], {i, 2, Length[data]}]]

Comments
Post a Comment