Variance

##sum of squares

#########Variance: A Worked Example##########
ozone<- span=""> read.table('gardens.txt',header=T)
attach(ozone)
ozone
##    gardenA gardenB gardenC
## 1        3       5       3
## 2        4       5       3
## 3        4       6       2
## 4        3       7       1
## 5        2       4      10
## 6        3       4       4
## 7        1       3       3
## 8        3       5      11
## 9        5       6       3
## 10       2       5      10
mean(gardenA)
## [1] 3
mean(ozone$gardenA)
## [1] 3
mean(ozone[,1])
## [1] 3
gardenA-mean(gardenA)
##  [1]  0  1  1  0 -1  0 -2  0  2 -1
(gardenA-mean(gardenA))^2
##  [1] 0 1 1 0 1 0 4 0 4 1
sum((gardenA-mean(gardenA))^2) #sum of squares
## [1] 12
sum((gardenA-mean(gardenA))^2)/(length(gardenA-1)) #the sum of squares #degree of 9
## [1] 1.2
mean(gardenB)
## [1] 5
gardenB-mean(gardenB)
##  [1]  0  0  1  2 -1 -1 -2  0  1  0
(gardenB-mean(gardenB))^2
##  [1] 0 0 1 4 1 1 4 0 1 0
sum((gardenB-mean(gardenB))^2)
## [1] 12
sum((gardenB-mean(gardenB))^2)/(length(gardenB)-1)
## [1] 1.333333
mean(gardenC)
## [1] 5
gardenC-mean(gardenC)
##  [1] -2 -2 -3 -4  5 -1 -2  6 -2  5
(gardenC-mean(gardenC))^2
##  [1]  4  4  9 16 25  1  4 36  4 25
sum((gardenC-mean(gardenC))^2)
## [1] 128
sum((gardenC-mean(gardenC))^2)/(length(gardenC)-1)
## [1] 14.22222
#F test 
var(gardenC)/var(gardenB)
## [1] 10.66667
2*(1-pf(10.667,9,9)) #The F Distribution
## [1] 0.001624002
help(pf)
## starting httpd help server ... done
var.test(gardenB,gardenB) #built-in F test:
## 
##  F test to compare two variances
## 
## data:  gardenB and gardenB
## F = 1, num df = 9, denom df = 9, p-value = 1
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2483859 4.0259942
## sample estimates:
## ratio of variances 
##                  1
gardenC
##  [1]  3  3  2  1 10  4  3 11  3 10
#########Variance and Sample Size
plot(c(0,32),c(0,15),type='n',xlab='Simple size',ylab='Variance')
for(n in seq(3,31,2)){
  for(i in 1:30){
    x<- span=""> rnorm(n,mean=10,sd=2)
    points(n,var(x))
  }
}


#########A Measure of Unreliability
sqrt(var(gardenA)/10)
## [1] 0.3651484
sqrt(var(gardenB)/10)
## [1] 0.3651484
sqrt(var(gardenC)/10)
## [1] 1.19257
########Confidence Intervals
#The first argument in qt() is the probability and the second is the #degrees of freedom
#left 95% confidence interval  value of Student’s t associated with α  
qt(.025,9)#Sample size 10,d.f 10-1
## [1] -2.262157
qt(.025,df=9)
## [1] -2.262157
##right 95% confidence interval  value of Student’s t associated with α 
qt(.975,9)
## [1] 2.262157
#values of t for 99% are bigger than these (0.005 in each tail):
qt(.995,9)
## [1] 3.249836
##99.5% confidence are bigger still (0.0025 in each tail):
qt(.9975,9)
## [1] 3.689662
#Bootstrap

data<- span=""> read.csv('skewdata.txt')
attach(data)
names(data)
## [1] "values"
#This draws the plain graph
plot(c(0,30),c(0,60),type='n',xlab='Sample size',ylab='Confidence interval')
for(k in seq(5,30,3)){
  a<- span=""> numeric(10000)
  for(i in 1:10000){
    a[i]<- span=""> mean(sample(values,k,replace=T))
  }
points(c(k,k),quantile(a,c(.025,.975)),type='b',pch=21,bg='red')
}


#At n= 30, the bootstrapped CI based on 10 00 simulations
quantile(a,c(.025,.975))
##     2.5%    97.5% 
## 24.90652 37.90855
#bootstrapped intervals compared with the intervals calculated from the normal (blue solid line
xv <- span=""> seq(5,30,0.1)
yv <- span=""> mean(values)+1.96*sqrt(var(values)/xv)
plot(c(0,30),c(0,60),type='n',xlab='Sample size',ylab='Confidence interval')
for(k in seq(5,30,3)){
  a<- span=""> numeric(10000)
  for(i in 1:10000){
    a[i]<- span=""> mean(sample(values,k,replace=T))
  }
  points(c(k,k),quantile(a,c(.025,.975)),type='b',pch=21,bg='red')
}
lines(xv,yv,col="blue")
yv <- span=""> mean(values)-1.96*sqrt(var(values)/xv)
lines(xv,yv,col="blue")
#Student’s t distribution (green dotted line):
yv <- span=""> mean(values)-qt(.975,xv-1)*sqrt(var(values)/xv)
lines(xv,yv,lty=2,col="green")
yv <- span=""> mean(values)+qt(.975,xv-1)*sqrt(var(values)/xv)
lines(xv,yv,lty=2,col="green")


#########################
points<- span=""> c(28,26,10,27,20,38,23,28,25,2)
x.m<- span=""> function(x) {
  #sprintf('x-mean=%f ',x-mean(x));
  x-mean(x)
}
x.m.t<- span=""> function(x) {
  #x.m() 
  #sprintf('x-mean(x)/length(x)=%f ',x-mean(x)/length(x))  
  ((x-mean(x))^2)
}
x_mean<- span=""> x.m(points)
x_mean_sqrt<- span=""> x.m.t(points)
x_mean
##  [1]   5.3   3.3 -12.7   4.3  -2.7  15.3   0.3   5.3   2.3 -20.7
x_mean_sqrt
##  [1]  28.09  10.89 161.29  18.49   7.29 234.09   0.09  28.09   5.29 428.49
sum(x_mean)/length(points)
## [1] 7.105427e-16
sum(x_mean_sqrt)/length(points)
## [1] 92.21
sqrt(sum(x_mean_sqrt)/length(points))
## [1] 9.602604
mean(points)
## [1] 22.7
(mean(points))^2
## [1] 515.29
var(points)
## [1] 102.4556