Monday 22 April 2013

Using R to test for a single Normal variance

Say we have a random sample from a Normal distribution whose variance, SigmaSquared, is unknown.

We may want to test the hypothesis that the population variance has a specific value, Sigma0Squared. The sample variance is SSquared, and the sample size is n.

One-sided test
To perform a one-sided test for the variance, we can use the R function oneSidedTestForNormalVariance() to perform a one-sided test:
> oneSidedTestForNormalVariance <- function(Sigma0Squared, SSquared, n, alternative)
   {
         T <- (n - 1)*(SSquared)/Sigma0Squared
         df <- n - 1
         if (alternative == 'less') { PValue <- pchisq(T, df)      }
         else                                { PValue <- 1 - pchisq(T, df) } 
         print(paste("Test statistic=",T,"and p-value=",PValue))
   }

This calculates a p-value for a one-sided test. For example, say we have a sample of 10 data points: 455.6, 457.5, 454.6, 456.9, 455.7, 451.6, 454.5, 456.7, 454.5, 451.3. The sample variance is 4.38. We can test the null hypothesis that the population variance has value 1.2, compared to the alternative hypothesis that the variance is greater than 1.2:
> x <- c(455.6, 457.5, 454.6, 456.9, 455.7, 451.6, 454.5, 456.7, 454.5, 451.3)
> oneSidedTestForNormalVariance(1.2, var(x), length(x), 'greater')
[1] "Test statistic= 32.8241666666664 and p-value= 0.000143299324329993"

Here the P-value is about 0.0001. Thus, there is strong evidence that the variance is not 1.2.

Two-sided test
In some cases, we may want to perform a two-sided test, then we can use this function:
> twoSidedTestForNormalVariance <- function(Sigma0Squared, SSquared, n)
   {
         T <- (n - 1)*(SSquared)/Sigma0Squared
         df <- n - 1
         if (Sigma0Squared > SSquared)      { PValue <- 2 * pchisq(T, df)         }
         else                                                   { PValue <- 2 * (1 - pchisq(T, df)) }  
         print(paste("Test statistic=",T,"and p-value=",PValue))
   }

For example, say the sample size (n) is 40, the sample variance (SSquared) is 165, and we want to perform a two-sided test of the hypothesis that the population variance is 225:
> twoSidedTestForNormalVariance(225, 165, 40)
[1] "Test statistic= 28.6 and p-value= 0.220604657817947"

Here the data do not provide sufficient evidence to reject the null hypothesis that the population variance is 225.

No comments: