C:\Rmisc>C:\R-2.9.0\bin\Rterm --vanilla R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > pdf("philrain.pdf") > > philrain <- read.fwf("C:\\data\\philrain.dat", width=c(8,8,8,8), + col.names=c("year","ymonth","smonth","rain")) > > philrain <- data.frame(philrain) > attach(philrain) > > ymonth <- factor(ymonth) > > # Summary Statistics for Rainfall > mean(rain) [1] 367.6796 > median(rain) [1] 340 > var(rain) [1] 36776.90 > sd(rain) [1] 191.7731 > min(rain) [1] 19 > quantile(rain,c(0.25,0.5,0.75)) 25% 50% 75% 232.75 340.00 468.00 > max(rain) [1] 1582 > > # POPULATION Mean and Variance saved to variables > MU.rain <- mean(rain) > SIGMA2.rain <- (length(rain)-1)*var(rain)/length(rain) > > > # Summary Statistics for rainfall by Month of Year > tapply(rain,ymonth,mean) 1 2 3 4 5 6 7 8 336.2889 300.5333 348.2000 365.0889 405.9111 417.6000 384.3333 429.2000 9 10 11 12 362.3556 340.0889 351.9556 370.6000 > tapply(rain,ymonth,sd) 1 2 3 4 5 6 7 8 175.8682 130.1317 149.7016 181.7714 192.9197 199.0010 183.8897 286.9646 9 10 11 12 248.8752 177.4534 146.0495 155.4938 > > # Plots of Rainfall > # Histogram > hist(rain,breaks=seq(0,1600,100)) > # Boxplot by Month of Year > plot(ymonth,rain) > # Time series > plot(rain,type="o",xlab="Month", ylab="Rainfall (1/100 inches)",main="Philadelphia Monthly Rain -- 1825-1869") > abline(h=mean(rain)) > > mean <- numeric(10000) > var <- numeric(10000) > # 10000 Random Samples of size n=30 (Sampling without replacement) > for (i in 1:10000) { + ransamp <- sample(rain,30) + mean[i] <- mean(ransamp) + var[i] <- var(ransamp) + } > > # Summary statistics for sample Mean and Variance > mean(mean) [1] 367.7721 > var(mean) [1] 1162.386 > sd(mean) [1] 34.09378 > mean(var) [1] 36751.92 > var(var) [1] 235190435 > sd(var) [1] 15335.92 > > min(mean) [1] 257.4333 > max(mean) [1] 496.7667 > min(var) [1] 9099.264 > max(var) [1] 126039.4 > > > > > # Histogram of Sampling Distribution of Y-bar with Normal Super-Imposed > hist(mean,breaks=seq(200,550,10),xlab="ybar",main="Sampling Distribution of Y-bar") > lines(225:515,10000*dnorm(225:515,MU.rain,sqrt(SIGMA2.rain/30))) > > # Histogram of (n-1)S^2/SIGMA^2 with Chi-squared Super-Imposed > hist(29*var/SIGMA2.rain,breaks=seq(0,200,2),xlab="S2",main="Sampling Distribution of Scaled S2") > lines(seq(0,200,2),10000*dchisq(seq(0,200,2),29)) > > > > dev.off() null device 1 >