统计软件选讲

目录Chapter 1-2Recipe1.Vector向量定义、性质、选取子集、判断Recipe2.Factor因子Recipe3.Matrix矩阵定义、转化、基本运算Recipe4.List清单Recipe5.Data Frame数据框Chapter 3Recipe1.Sample随机数Recipe2.Central Limit Theorem中心极限定理Chapter 4Recipe1.Functions函数Recipe2.Loops循环Forloops、ifelse、while、stretch、powRecipe3.防错&除错Max.power、max.power_bookRecipe4.Examples例子Basic_funs、is.symChapter 5Recipe1.读取和整理Chapter 6Recipe1.Linear Regression Model线性回归模型散点图、回归分析、作图、分析图，多元线性回归、lmcurve Recipe2.Analysis of Variance Model方差分析模型盒子图、方差分析Recipe3.Generalized Linear Models广义线性模型Binomial Model(Logistic Regression)Chapter 7Recipe1.Linear Discriminant Analysis线性判别分析Recipe2.Cluster Analysis聚类分析分层聚类分析、非分层聚类分析Recipe3.Library程序馆Chapter 8Recipe1.Multiframe Graphic多重图框Recipe2.增加直线一、R软件教程Chapter 1-2Recipe1.Vector向量rm(list = ls())x <- c(1, 2, 3, 4) ##c(向量内容，若是字符需要加"")y <- 1:100 ##连续型z <- seq(1, 100, by = 8) ##等差数列seq(首项, 最大值, by = 公差)w <- rep(1:4, each = 2, times = 2)##重复性rep(内容, c(各项重复次数) or each = n全部都重复n次, times = m整体重复m次)#####3倍向量#####r <- NULLfor (i in 1:50){r <- c(r,3*i)}#################性质length(x) ##求长度mean(y) ##求均值var(z) ##求方差#选取子集x[1] ##选取第一个分量y[2:4] ##选取第2到4个分量z[c(1, 5, 1)] ##选取第1，5，1个分量w[-1] ##去除第1个分量r[r < 25] ##选取小于25的分量#判断Is.vector(x) ##是否为向量Is.numeric(x) ##是否为数字构成的向量Recipe2.Factor因子rm(list = ls())medicine <- c(rep("Jack", 3), rep("A", 4), rep("B", 3))is.factor(medicine) ##检验is.**()medicine <- as.factor(medicine) ##转化as.**() 因子levels(medicine) ##水平Recipe3.Matrix矩阵rm(list = ls())x1 <- matrix(1:9, 3, ) ##（向量，行数，列数），向量数=行数*列数，默认按列排列x2 <- matrix(2:10, , 3, byrow=TRUE) ##byrow=TRUE按行排列diag(5) ##5*5单位阵dim(x1) ##维数colnames(x1) <- c("c1", "c2", "c3") ##列名rownames(x1) <- c("r1", "r2", "r3") ##行名x2[c(1, 3), c(2, 3)] ##调动子矩阵#转化t1 <- as.vector(x1) ##将x1转变为向量t2 <- as.matrix(t1) ##将t1转变为矩阵，默认列向量t(x1) ##转置t(t1) ##将t1转为矩阵，默认行向量q <- solve( matrix( c(1, 5, 6, 3), 2, ) ) ##求逆矩阵dim(q) ##求行数列数#基本运算x1 + x2 ##矩阵加法x1 * x2 ##对应元素相乘（错误做法，注意）x1 %*% x2 ##矩阵乘法v1 <- 1:3x1 * v1 ##第n行*第n元素（错误做法，注意）x1 %*% as.matrix(v1) ##向量乘矩阵要先转化A <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 10), 3, , byrow = TRUE)e <- eigen(A) ##求特征值及特征向量e$values ##调用特征值svd(A) ##奇异值分解B <- matrix(1:6, 2, )C <- matrix(5:10, , 2)rbind(A, B) ##按行合并cbind(A, C) ##按列合并apply(A, 1, sum) ##按行加和apply(A, 2, sum) ##按列加和sd(标准差）medina####定义3维矩阵####a <- 1:24dim(a) <- c(2, 3, 4)#################Recipe4.List清单rm(list = ls())a <- 1b <- 5:8c <- matrix(1:10, 2, )d <- c("Jack", "Lucy", "Mark", "L.C.")x <- list(sex = a, age = b, birth = c, name = d) ##定义清单各名字x[[1]] ##调用第1对象x$age ##调用此名称对象x[["birth"]] ##调用此名称对象Recipe5.Data Frame数据框rm(list = ls())grade <- read.csv("Book1.csv") ##读取数据part <- names(grade) ##名字grade <- as.matrix(grade)grade <- as.factor(grade)gradepart <- rep(part, each=16)parttemp <- matrix(c(grade, part), nc=2) ##结合grade与part，nc列数temptemp <- temp[!is.na(temp[,1])] ##除去NA项temptemp <- matrix(c(grade, part), nc=2)tempbeauties <- data.frame(temp) ##数据框beautiescolnames(beauties) <- c("Grade", "Name") ##改名字attach(beauties) ##通过名字寻访数据###########################################rm(list=ls())grade <- read.csv("Book1.csv") ##读取数据gradegrade1 <- read.csv("Book1.csv", header=FALSE) ##将名字默认为数据grade1attributes(grade) ##总结数据grade[,1]min(grade[,1])grade[2,]max(grade[2,])attach(grade) ##方便调用g_matrix <- as.matrix(grade) ##数据框转矩阵g_matrixis.matrix(g_matrix)g_frame <- data.frame(g_matrix) ##矩阵转数据框g_frameis.matrix(g_frame)Recipe1.Sample随机数rm(list = ls())set.seed(20140919) ##设置随机种子，下次再启用数据不变v1 <- sample(1:100, 9, replace = TRUE) ##随机向量，有放回v2 <- sample(1:100, 9, replace = FALSE)Recipe2.Central Limit Theorem中心极限定理rm(list = ls())set.seed(12345)u = runif(10000) ##均匀分布summary(u) ##总结var(u) ##方差hist(u, main = "U(0, 1)") ##直方图u = matrix(u, nc = 10)m = apply(u, 1, mean) ##计算u内每一行的平均数summary(m)var(m)hist(m, main = "样本平均(n=10)的抽样分布")Recipe1.Functions函数rm(list = ls())fix(ifelse) ##找出函数，修改后源文件不变，但再次fix版改变edit(max.power) ##找出函数，修改后皆不变Recipe2.Loops循环######forloops######求和rm(list = ls())VC_loops <- function(a){##a must be a vector!n <- length (a)s <- 0for(i in 1:n){ s <- s + a[i] }s ##or return (s)}————RUN————rm(list = ls())source("forloops.r")x <- c(1, 2, 3, rep(3, 10), sample(10))VC_loops(x)######ifelse######求数值为0和21的数的总数rm(list = ls())VC_ifelse1 <- function(vec){n <- length(vec)s <- 0for(i in 1:n){if(vec[i]==0){ s <- s + 1 }else if(vec[i]==21){ s <- s + 1 }}s}VC_ifelse2 <- function(vec){n <- length(vec)s <- 0for(i in 1:n){if(vec[i]==0 || vec[i]==21)##||or &&and ==is !=not{ s <- s + 1 }else{ s <- s }}return (s)}————RUN————source("ifelse.r")VCvec <- c(0, 56, 21, 21, 21, 0, 59, 98, 356, 847, 21, 0)VC_ifelse1(VCvec)VC_ifelse2(VCvec)######while######rm(list = ls())VC_while <- function(x, y){n <- 0while (x < y){n <- n + 1x <- x^n}n}————RUN————source("while.r")x <- 10y <- 1001VC_while(x, y)######stretch######rm(list = ls())VC_stretch <- function(vec, num){output <- NULLnn <- length(vec)for(i in 1:nn){ output <-c(output, rep(vec[i], num)) } output}————RUN————source("stretch.r")VCvec2 <- c("VC", "JC", "Cara", "Cuckoo")VC_stretch(VCvec2, 10)######pow######rm(list = ls())VC_pow <- function(x, y){if(x > y){ return ("Error! x is bigger than y!") } else{n <- 0for(i in 1:100){if(x^(n+1) < y){ n <- n + 1 }}x^n}}————RUN————source("pow.r")x <- 10y <- 1001z <- 101VC_pow(x, y)VC_pow(y, z)Recipe3.防错&除错rm(list = ls())######max.power#######max.power <- function(x, y){if(x < y & x > 0){n <- 1midx <- xwhile(x < y){n <- n + 1x <- midx ^ n}return(n - 1)}Else{ stop("x is bigger than y or something else!") ##防错} }debug(max.power) ##除错max.power(2,9) ##运行后，按步骤运行undebug(max.power) ##退出除错######max.power_book######max.power_book <- function(x, y){if(x < y & x > 0)if(x > 1){i <- 1total <- xwhile(total*x < y){i <- i+1total <- total*x}i}else"Infinity"else NULL}Recipe4.Examples例子######Basic_funs######rm(list = ls()) ##Definition of functions in R! VC_sum <- function(a, b){ a+b } ##定义VC_sumVC_odds <- function(a, b) ##定义VC_odds {if(a-b !=0){ return((a+b)/(a-b)) }else{ return("Error!") }}————RUN————rm(list = ls())source("Basic_funs.r") ##读取文件a <- 2b <- 2VC_sum(a, b)VC_odds(a, b)######is.sym######rm(list = ls())is.sym <- function(x){if(!is.matrix(x)){ stop("The input must be a matrix!") }n <- dim(x)[1]m <- dim(x)[2]if(n != m){ stop("Matrix is not squared!") }for(i in 1:n){for(j in 1:n){if(x[i, j] != x[j, i]){ stop("Matrix is not synetric!") } }}return("The input is a symetric matrix!")}————RUN————rm(list = ls())source("is.sym.r")x <- 1:9is.sym(x)x <- matrix(1:9, 3, )is.sym(x)x <- matrix(1:12, 3, )is.sym(x)x <- matrix(rep(1, 9), 3, )is.sym(x)Recipe1.读取和整理rm(list = ls())x <- read.csv("replace.csv") ##读取数据y <- replace(x, is.na(x), 0) ##替换数据z <- replace(y, y < 6, NA)Recipe1.Linear Regression Model线性回归模型rm(list = ls())data(faithful) ##调动R中原有数据attach(faithful) ##分量数据化plot(eruptions, waiting) ##散点图faithful.lm <- lm(waiting ~ eruptions) ##回归分析summary(faithful.lm) ##总结curve(33.47 + x * 10.73, 1.5, 5.0, add = TRUE, lwd = 3, col = 4)##作图(方程，左区间，右区间，add:if TRUE 加入已有的散点图，,线粗，颜色)par.lm <- par(mfrow = c(2,2))plot(faithful.lm) ##分析图par(par.lm)######多元线性回归######rm(list = ls())data(stackloss)attach(stackloss)y <- stackloss&stack.lossx1 <- Air.Flowx2 <- Water.Tempx3 <- Acid.Conc.a <- lm(y ~ x1 + x2 + x3)summary(a)b <- lm(y ~ x1 + x2 + x3, subset = -21)summary(b)######lmcurve######rm(list = ls())x <- rnorm(100)y <- rnorm(100)a <- lm(y ~ x)res <- a$coefficientsplot(x, y)curve(res[1] + res[2]*x,-3, 3)rm(list = ls())for(i in 1:3){x <- rnorm(10)y <- rnorm(10)a <- lm(y ~ x)res <- a$coefficientsplot(x, y)curve(res[1] + res[2]*x,-3, 3, lwd = 2, col = 3, add = T) savePlot(paste(i), type = "jpg")}Recipe2.Analysis of Variance Model方差分析模型rm(list = ls())data(PlantGrowth)attach(PlantGrowth)plot(weight~group, main = "Plant weight", xlab = "Fri")##盒子图黑线为中位数，框上下线分别为75%、25%分位数，上下线分别为最大小值————————————————x_weight <- weightx_group <- groupx_group <- as.vector(x_group)x_group[which(x_group == "ctrl")] <- 0x_group[which(x_group == "trt1")] <- 1x_group[which(x_group == "trt2")] <- 2windows()plot(x_group, x_weight, main = "Plant weight", xlab = "Fri")##与上图作对比，如果x轴为数字，一定要把它变为因子————————————————windows()x-group <- as.factor(x_group)plot(x_group, x_weight, main = "Plant weight", xlab = "Fri")a <- aov(weight ~ group) ##方差分析summary(a)Recipe3.Generalized Linear Models广义线性模型####Binomial Model(Logistic Regression)####rm(list = ls())neuro <- read.csv("neuro.csv")attach(neuro)y <- cbind(Cured, NotCured)a <- glm(y ~ Sex + Treat, family = binomial) ##广义线性模型anova(a)Recipe1.Linear Discriminant Analysis线性判别分析rm(list = ls())library(MASS) ##程序馆data(iris)attach(iris)a <- lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width)##线性判别分析aRecipe2.Cluster Analysis聚类分析rm(list = ls())set.seed(20141112)x <- rnorm(100)a <- hclust(dist(x)) ##分层聚类分析plot(a)kmeans(x, 3) ##非分层聚类分析Recipe3.Library程序馆rm(list = ls())##以lasso2为例##download lasso2 from Beijing2 first!install.packages("lasso2")library(lasso2)data(Iowa)l1c.I <- l1ce(Yiels ~ ., Iowa, bound = 10, absolute.t = TRUE)l1c.IRecipe1.Multiframe Graphic多重图框rm(list = ls())data(faithful)attach(faithful)op <- par(mfrow = c(2, 2)) ##2行2列hist(waiting) ##直方图hist(eruptions)qqnorm(waiting)qqline(waiting) ##qqnorm与qqline在一个图里qqnorm(eruptions)qqline(eruptions)par(op) ##停止作图savePlot("hist_qq", type = "jpeg")Recipe2.增加直线rm(list = ls())data(faithful)attach(faithful)plot(eruptions, waiting, main = "Springs")abline(lsfit(eruptions, waiting), lwd=3, col=5) ##拟合直线abline(h=70, v=3, lwd=2, col=4) ##随便画的，h平行于x轴，v平行于y轴。