options(repos=c(CRAN="http://mirrors.ustc.edu.cn/CRAN/")) #mathematical statistics examples with R #----------------------------------- #summary statistics #----------------------------------- x<-rnorm(20,0,2) summary(x) # moments library(moments) moment(x,order= 2, central=T) #sample mean, median and variance mean(x); median(x); var(x) #skewness and kurtosis skewness(x); kurtosis(x) #empirical cdf x<-rbinom(10,8,0.2) Fn<-ecdf(x) plot(Fn) #-------------------------------- # Estimates #-------------------------------- #MLE #first, get sample xvec <- c(2,5,3,7,-3,-2,0) #second, write a function which is negative of the log likelihood fn <- function(theta) { sum ( 0.5*(xvec - theta[1])^2/theta[2] + 0.5* log(theta[2]) ) } #where there are two parameters: theta[1] and theta[2]. #find the mle nlm(fn, theta <- c(0,1), hessian=TRUE) #--------------------------------- # hypothesis testing #--------------------------------- # one-sample t test x<-c(-1.892,2.791,3.280,0.911,-0.682,-0.646,1.240,0.557, 0.323, -4.138) t.test(x) # two-sample t test x<-c(-1.892,2.791,3.280,0.911,-0.682,-0.646,1.240,0.557, 0.323, -4.138) y<-c(-1.457,3.172,2.167,5.480,-0.021,-1.847,-5.712,1.588, 1.055,-0.319) t.test(x,y,var.equal=TRUE) #or you can do var equality test var.test(x,y) # chi-square test heartatk<-read.table("heartatk4R.txt",header=TRUE) dim(heartatk) attach(heartatk) chisq.test(table(SEX)) table(SEX,DIED) chisq.test(table(SEX,DIED)) #kolmogorov-smirnov test x<-c(-1.892,2.791,3.280,0.911,-0.682,-0.646,1.240,0.557, 0.323, -4.138) y<-c(-1.457,3.172,2.167,5.480,-0.021,-1.847,-5.712,1.588, 1.055,-0.319) ks.test(x,pnorm) ks.test(x,y) wilcox.test(x,y) z<-x-y np<-sum(z>0) binom.test(np,length(x)) #test of normality shapiro.test(x) qqnorm(x) library(nortest) library(help=nortest) # lienar regression library(faraway) data(savings) g <- lm(sr~pop15 + pop75 + dpi + ddpi, data=savings) summary(g) #testing H0: b1=b2=b3=b4=0 rssr<-sum((savings$sr-mean(savings$sr))^2) rssf<-sum(g$res^2) n<-nrow(savings) p<-ncol(savings)-1 dfr<-n-1 dff<-n-p-1 ft<-((rssr-rssf)/p)/(rssf/(n-p-1)) 1-pf(ft,p,n-p-1) #use anova gr<-lm(sr~1, data=savings) anova(gr,g) #Testing just one predictor # H0: b1=0 g2 <- lm(sr~pop75 + dpi + ddpi, data=savings) rssr2<-sum(g2$res^2) ft2<-(rssr2-rssf)/(rssf/(n-p-1)) 1-pf(ft2,1,n-p-1) #use anova anova(g2,g) #Testing a subspace #H0 : bpop15 = bpop75 gr <- lm(sr~I(pop15+pop75)+dpi+ddpi,savings) anova(gr,g) #H0 : bddpi=1 gr <- lm(sr~pop15+pop75+dpi+offset(ddpi),savings) anova(gr,g)