R语言三分位 r语言求四分位极差

data.entry(mtcar) # 编辑
 edit(mtcar) # 编辑
 fix(mtcar) # 列出结构
 attach(mtcar)
 detach(mtcar)
 table(mtcar) 
 barplot(table(Cry))
 mean(mtcars$mpg,trim=0.1)  # 截外平均
 mean(mtcars$mpg）  # 平均
 tapply(mtcars$mpg,mtcars$cyl,mean)  # 按组平均
 mean(mtcars$mpg[mtcars$cyl==4])  # 条件平均
 IQR(mpg) # 四分位数的极差
  quantile(mpg) ； fivenum(mpg)  分位数
 summary(mtcar)
 plot(cyl,mpg)
 plot(hp,mpg,pch=cyl)
 legend(250,30,pch=c(4,6,8),legend=c("4 cylinders","6 cylinders","8 cylinders"))
  z<-lm(cyl~mpg)
  cor(cyl,mpg) # 相关系数 注意有至少三种相关系数
  lm.resids <- resid(lm.res) # 残差
  plot(lm.resids)
  hist(lm.resids)
  qqnorm(lm.resids)
  mode(x) # 数据类型
  Inf  -Inf  # 无限
   x<-5/0
   ls( )
  ls.str( )
  ls.str(pat="M", max.level=-1)
  rm(list=ls( ))
  rm(list=ls(pat="^m")) # 删除内存中以m开头的对象
  z <- seq(1,5,by=0.5) 
  z <- seq(1,10,length=11)
  z <- rep(1:3, times = 4, each = 2)
  z <- rep(2:5,rep(2,4))
  rep(2:6,rep(1:5))
  rep(1:3, times = 4, each = 2)
  rep(1:4, c(2,1,2,1))
  sequence(3:5)
  sequence(c(10,5))
  labs <- paste(c("X","Y"), 1:10, sep="")
  x <- c(10.4, 5.6, 3.1, 6.4, 21.7)
   temp <- x > 13
   a <- c("green", "blue", "green", "yellow")
   a <- factor(a)
   ff <- factor(c("A", "B", "C"), labels=c(1,2,3))
   
   b <- c(1,2,3,1)
  b <- factor(b)
 levels(b) <- c("low", "middle", "high")
 f <- factor(c(2, 4), levels=2:5)
  gl(3, 5)
   x <- c(10.4, 5.6, 3.1, 6.4, 21.7) # 向量运算
   y <- c(x,0,x)
    v <- 2*x + y + 1  # 向量运算，循环法则，要整数陪
    which.max（X） # 返回向量中最大值的下标
    median（X）  # 中位数
    rev(x)
    rank(x)
    var(x,y)
    cov(x,y)
    cor(x,y)
    x[c(1，4)]
    x[-(2:5)]
    fruit <- c(5, 10, 1, 20)
     names(fruit) <- c("orange", "banana", "apple", "peach")
x[x>10]
y = runif(100,min=0,max=1)
y <- x[!is.na(x)]
dim(x)
A <- array(1:8, dim = c(2, 2, 2))
dimnames(A) <- list(c("a", "b"), c("c", "d"), c("e", "f"))
colnames(A)
rownames(A)
matrix()
diag(x)
diag(2.5, nr = 3, nc = 5)
X <- matrix(1:4, 2) # 等价于X <- matrix(1:4, 2, 2)
X <- matrix(1:4, 2, 4, byrow=TRUE) #按行填充
x[2,2]
x[2,]
x[,c(2,3),drop=FALSE]
x[,-2]
 x[,3] <- NA
     x[is.na(x)] <- 1
t(x)
rbind(x1,x2)
cbind(x1,x2)
m2*m2  # 矩阵乘集
rbind(m1, m2) %*% cbind(m1, m2)  #%*% 矩阵的代数乘积
 apply(X, MARGIN, FUN)  #MARGIN=1表示按列计算, MARGIN=2表示按行计算, MARGIN=c(1,2)表示按行列计算(在至少3维的数组中使用)
 class(F), mode(F)
 scale(m, center=T, scale=T)
 z.df=data.frame(INDEX = y, VALUE = x)
 foo <- read.table(file = "c:/data/foo.txt", header = T)
 sweep(x,1,STATS=row.med,FUN="-") #可以设置减去一个值
 pairs(Puromycin, panel = panel.smooth)
 xtabs(~state + conc, data = Puromycin)  #交叉因子产生的列联表
 subset(Puromycin, state == "treated" & rate > 160)
 subset(Puromycin, conc > mean(conc))
 Puromycin$iconc <- 1/Puromycin$conc
 Puromycin$iconc <- with(Puromycin, 1/conc)
 Puromycin <- transform(Puromycin, iconc = 1/conc,sqrtconc = sqrt(conc)) 
 head(Puromycin)
 L1 <- list(1:6, matrix(1:4, nrow = 2))
 L2 <- list(x = 1:6, y = matrix(1:4, nrow = 2))
 ts(1:47, frequency = 12, start = c(1959, 2)) #时间序列
 setwd("C:/data")
  write.table(d, file = "c:/data/foo.txt",row.names = F, quote = F)
  write.csv(d, file = "c:/data/foo.csv",row.names = F, quote = F)  
  save(list =ls(all=TRUE), file=".RData")
  HousePrice <- read.table(file="houses.dat")
  HousePrice <- read.csv("houses.csv", header=TRUE)
   mydata <- scan("data.dat", what = list(Sex="", Weight=0, Height=0)) #读取文档，分配类型
mydata <- read.fwf("data.txt", widths=c(1, 4, 3),col.names=c("X","Y","Z")) #读固定宽度文件
mydata <- read.delim("clipboard")
library(RODBC)
          z <- odbcConnectExcel("c:/data/body.xls")
          foo <- sqlFetch(z, "Sheet1")
          close(z)
data( )
data(package="pkname")
data(dataname, package="pkname")
mtcars2 <- data.frame(mtcars[,c(1,4)])
load("c:/data/myR.Rdata")
demo(graphics）
plot(x, y)
pie(x) 饼图
 boxplot(x) 盒形图(\box-and-whiskers")
 stripchart(x)  把x的值画在一条线段上， 样本量较小时可作为盒形图的替代
 coplot(x~y | z) 关于z的每个数值(或数值区间)绘制x与y的二元图
 interaction.plot (f1, f2, y) 如果f1和f2是因子， 作y的均值图， 以f1的不同值作为x轴， 而f2的不
 同值对应不同曲线； 可以用选项fun指定y的其他的统计量(缺省计算
 均值， fun=mean)
 matplot(x,y) 二元图， 其中x的第一列对应y的第一列， x的第二列对应y的第二列，
 依次类推.
 dotchart(x) 如果x是数据框， 作Cleveland点图(逐行逐列累加图)
 fourfoldplot(x) 用四个四分之一圆显示2times2列联表情况(x必须是dim=c(2, 2,
 k)的数组， 或者是dim=c(2, 2)的矩阵， 如果k  1)
 assocplot(x) Cohen{Friendly图， 显示在二维列联表中行、 列变量偏离独立性的
 程度
 mosaicplot(x) 列联表的对数线性回归残差的马赛克图
 pairs(x) 如果x是矩阵或是数据框， 作x的各列之间的二元图
 plot.ts(x) 如果x是类"ts"的对象， 作x的时间序列曲线， x可以是多元的， 但是
 序列必须有相同的频率和时间
 ts.plot(x) 同上， 但如果x是多元的， 序列可有不同的时间但须有相同的频率
 hist(x) x的频率直方图
 barplot(x) x的值的条形图
 qqnorm(x) 正态分位数－分位数图
 qqplot(x, y) y对x的分位数－分位数图
 contour(x, y, z) 等高线图(画曲线时用内插补充空白的值)， x和y必须为向量， z必须为
 矩阵， 使得dim(z)=c(length(x), length(y)) (x和y可以省略)
 filled.contour (x, y, z) 同上， 等高线之间的区域是彩色的， 并且绘制彩色对应的值的图例
 image(x, y, z) 同上， 但是实际数据大小用不同色彩表示
 persp(x, y, z) 同上， 但为透视图
 stars(x) 如果x是矩阵或者数据框， 用星形和线段画出
 symbols(x, y, ...) 在由x和y给定坐标画符号(圆， 正方形， 长方形， 星， 温度计式或者盒
 形图)， 符号的类型、 大小、 颜色等由另外的变量指定
 termplot(mod.obj) 回归模型(mod.obj)的(偏)影响图
 points(x, y) 添加点(可以使用选项type=)
 lines(x, y) 同上， 但是添加线
 text(x, y, labels,
 ...)
 在(x,y)处 添 加 用labels指 定 的 文 字； 典 型 的 用 法 是: plot(x, y,
 type="n"); text(x, y, names)
 mtext(text, side=3,
 line=0, ...)
 在边空添加用text指定的文字， 用side指定添加到哪一边(参照下面的axis(
 ))； line指定添加的文字距离绘图区域的行数
 segments(x0, y0, x1,
 y1)
 从(x0,y0)各点到(x1,y1)各点画线段


 arrows(x0, y0,
 x1, y1, angle= 30,
 code=2)
 同上, 但加画箭头. 如果code=2, 则在各(x0,y0)处画箭头; 如果code=1, 则在
 各(x1,y1)处画箭头; 如果code=3, 则在两端都画箭头angle控制箭头轴到箭头
 边的角度.
 abline(a,b) 绘制斜率为b和截距为a的直线
 abline(h=y) 在纵坐标y处画水平线
 abline(v=x) 在横坐标x处画垂直线
 abline(lm.obj) 画由lm.obj确定的回归线
 rect(x1, y1, x2, y2) 绘制长方形， (x1, y1)为左下角， (x2,y2)为右上角
 polygon(x, y) 绘制连接各x,y坐标确定的点的多边形
 legend(x, y, legend) 在点(x,y)处添加图例， 说明内容由legend给定
 title( ) 添加标题， 也可添加一个副标题
 axis(side, vect) 画坐标轴. side=1时画在下边; side=2时画在左边; side=3时画在上边;
 side=4时画在右边. 可选参数at指定画刻度线的位置坐标
 box( ) 在当前的图上加上边框
 rug(x) 在x-轴上用短线画出x数据的位置
 locator(n, type="n",
 ...)
 在用户用鼠标在图上点击n次后返回n次点击的坐标(x; y)； 并可以在点击处绘
 制符号(type="p"时)或连线(type="l"时)， 缺省情况下不画符号或连线


 text(x, y,expres\-sion(...))
 par(bg="yellow")
 head(Puromycin)


 plot(rate ~ conc, data = PuroA)
 with(PuroA, plot(conc, rate))
 plot(PuroA$rate, PuroA$conc)
  plot(u ~ 1, pch = u, col = u, cex = 3)
  
  plot(rate~conc,data=P,pch=2,col=4,cex=2.5,xlim=c(0,1.2),ylim=c(40,210), ylab="comm",xlab="rate",cex.lab=2)
  plot(rate ~ conc, data = P)
 > smooth1 <- with(P, lowess(rate ~ conc, f = 0.9))
 > lines(smooth1, col = "red")
 m2 <- lm(rate ~ conc + I(conc^2), data = P)
 lines(fitted(m2) ~ conc, data = PuroA, col = "blue") #增加线
 abline(lm(rate ~ conc, data = PuroA))
 mysymb <- c(1, 2)[Puromycin$state]
 plot(rate ~ conc, data = Puromycin, col = mysymb,
 pch = mysymb)
 legend(x = 0.6, y = 100,
 legend = c("treated", "untreated"),
 col = c(1, 2), pch = c(1, 2), lty = c(1, 3))
 par(mfrow = c(1, 2))
 dev.off( ) # 关闭原图
      layout(matrix(1:m*n, m, n)
x <- pmin(3, pmax(-3, stats::rnorm(50)))
 > y <- pmin(3, pmax(-3, stats::rnorm(50)))
 > xhist <- hist(x, breaks=seq(-3,3,0.5), plot=FALSE)
 > yhist <- hist(y, breaks=seq(-3,3,0.5), plot=FALSE)
 > top <- max(c(xhist$counts, yhist$counts))
 > xrange <- c(-3,3); yrange <- c(-3,3)
 > layout(matrix(c(2,0,1,3), 2, 2, byrow=TRUE),
 c(3,1), c(1,3), TRUE)
 > layout.show(3) # 给出作图窗口及编号
 > par(mar=c(3,3,1,1)) # 设定边界空行数
 > plot(x, y, xlim=xrange, ylim=yrange, xlab="", ylab="")
 > par(mar=c(0,3,1,1))
 > barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0)
 > par(mar=c(3,0,1,1))
 > barplot(yhist$counts, axes=FALSE, xlim=c(0, top),
 space=0, horiz=TRUE)
  ifelse(x >= 0, sqrt(x), NA)
  if (x >= 0) sqrt(x) else NA
  for(i in 4:10) print(1:i)
  while (i<5){print(1:i);i=i+1}
  for (i in 1:length(x)){
 if (x[i] == b)
 y[i] <- 0
 else
 y[i] <- 1
 }
 y[x == b] <- 0
  myfun <- function(S, F) {  # 编写函数
 data <- read.table(F)
 plot(data$V1, data$V2, type="l")
 title(S)
 }
 sapply(species, myfun, file)  #批量应用函数
 sample(c("成功", "失败"), 10, replace=T, prob=c(0.9,0.1)) # prob 出现的概率
 1/prod(52:49)  prod return all number's product
 1/choose(52,4) 从52个里面没有次序的选出4个 


 qnorm(0.025) 注意q, p, d, r 加载在不同分布前面的意思
 1 - pchisq(3.84, 1)
 2*pt(-2.43, df = 13)


 hist(x, breaks = "Sturges", freq = NULL, probability = !freq,
 hist( )的调用格式
 col = NULL,
 main = paste("Histogram of" , xname),
 xlim = range(breaks), ylim = NULL,
 xlab = xname, ylab,
 axes = TRUE, nclass = NULL)


 density(x, bw = "nrd0",
 density( )的调用格式
 kernel = c("gaussian", "epanechnikov", "rectangular",
 "triangular", "biweight", "cosine", "optcosine"),
 n = 512, from, to)
  fivenum(mtcars)
  quantile(mtcars$cyl)
  quantile(fpossum$totlngth ,prob=c(0.25,0.5,0.75))
  var(fpossum$totlngth)
  skewness(fpossum$totlngth)
  kurtosis(mtcars)
  plot(cars$speed~cars$dist,xlab="speed",ylab="dist")
  lines(lowess(cars$speed,cars$dist),lwd=3)
  rug(side=2, jitter(cars$dist, 20))
  contour(z, col = "red", drawlabels = FALSE,
 main = "Density estimation: contour plot")
  par(mfrow=c(1,2))
 > plot(brain~body,data=Animals)
 > plot(log(brain)~log(body),data=Animals)   # Log 变换
 >plot(d) # 或者 pairs(d) 多维数据
 matplot(d, type = 'l', ylab = "", main = "Matplot")
 aggregate(state.x77, list(Region = state.region), mean)
 aggregate(state.x77,list(Region=state.region,Cold=state.x77[,'Frost']>130),mean)
 options(digits=3)
 sd(state.x77)
 coplot(y ~ x | a)
 boxplot(length~species,data=cuckoos，xlab="length of egg",horizontal=TRUE))  #框须图
 stripchart(cuckoos$length~cuckoos$species, method ="jitter" ) # 条形图
 table( ), xtabs( ) 或ftable( )
 margin.table(Eye.Hair,1)
 prop.table
 a <- as.table(apply(HairEyeColor,c(1,2),sum))
 barplot(a, legend.text = attr(a, "dimnames")$Hair)
 barplot(a, beside = TRUE,
 legend.text = attr(a, "dimnames")$Hair)
 x <- 1:10
 attr(x,"dim") <- c(2, 5)
 optim( )或者nlm( )来求似然函数的极大值
 f <- function(P)(P^517)*(1-P)^483
 optimize(f,c(0,1),maximum = TRUE)  


  function(x,n,sigma,alpha,u0=0,alternative="two.sided"){
     options(digits=4)
     result<-list( )
     mean<-mean(x)
     z<-(mean-u0)/(sigma/sqrt(n))
     p<-pnorm(z,lower.tail=FALSE)
     result$mean<-mean
     result$z<-z
     result$p.value<-p
     if(alternative=="two.sided"){
         p<-2*p
         result$p.value<-p
     }
     else if (alternative == "greater"|alternative =="less" ){
         result$p.value<-p
     }
     else return("your input is wrong")
     result$conf.int<- c(
         mean-sigma*qnorm(1-alpha/2,mean=0, sd=1,
                          lower.tail = TRUE)/sqrt(n),
         mean+sigma*qnorm(1-alpha/2,mean=0, sd=1,
                          lower.tail = TRUE)/sqrt(n))
     result
 }






 conf.int<-function(x,n,sigma,alpha){
 options(digits=4)
 mean<-mean(x)
 c(mean-sigma*qnorm(1-alpha/2,mean=0, sd=1,
 lower.tail = TRUE)/sqrt(n),
 mean+sigma*qnorm(1-alpha/2,mean=0, sd=1,
 lower.tail = TRUE)/sqrt(n))
 }


 t.test(x)$conf.int


 two.sample.ci<-function(x,y,conf.level=0.95, sigma1,sigma2 ){
 options(digits=4)
 m= length(x); n = length(y)
 xbar=mean(x)-mean(y) alpha = 1 - conf.level
 zstar= qnorm(1-alpha/2)* (sigma1/m+sigma2/n)^(1/2)
 xbar +c(-zstar, +zstar)
 }
 var.test(x,y)
 prop.test(38,200,correct=TRUE)


 > like<-c(478, 246)
 > people<-c(1000, 750)
 > prop.test(like, people)


 source("shi.R")
  t.test(CaCo3, mu=20.7)
  var.test(x，y)
  t.test(x, y, paired=TRUE) #成对数据的检验


  onesamp( ) # DAAG包
  
  binom.test( )
  prop.test(35, 120, p=0.25, conf.level=0.975, correct=TRUE)
  
  sign.test<-function(x,m0,alpha=0.05,alter="two.sided"){
 +     p<-list()
 +     n<-length(x)
 +     sign<-as.numeric(x>=m0)
 +     s<-sum(sign)
 +     result<-binom.test(s,n,p=0.5,alternative=alter,conf.level=alpha)
 +     p$p.value=result$p.value
 +     p
 + }


      fisher.test(compare, alternative = "greater")  #列联表
 wilcox.test(diabetes,normal,exact = FALSE, correct=FALSE) # 非参数，是否相等
 mood.test(A,B) # 非参数，变化程度
 
 kruskal.test(x)
 ansari.test(worker.a,worker.b) # 尺度参数的检验
 fligner.test(x)
 
 X<-c(25.6, 22.2, 28.0, 29.8, 24.4, 30.0, 29.0, 27.5, 25.0, 27.7,
   23.0, 32.2, 28.8, 28.0, 31.5, 25.9, 20.6, 21.2, 22.0, 21.2)
   A<-factor(rep(1:5, each=4))
   miscellany<-data.frame(X, A)
   aov.mis<-aov(X~A, data=miscellany) 
   pairwise.t.test( ) # 两两检验
   pairwise.t.test(X, A, p.adjust.method="none")
   pairwise.t.test(X, A, p.adjust.method="holm")
   TukeyHSD(aov(X~A, sales))  #同时置信区间
   bartlett.test(X~A, data=sales) # 方差齐性检测
   library(car)
   levene.test(sales$X, sales$A) # 方差齐性检测
   oneway.test( )
   
   juice<-data.frame(
 X = c(0.05, 0.46, 0.12, 0.16, 0.84, 1.30, 0.08, 0.38, 0.4,
 0.10, 0.92, 1.57, 0.11, 0.43, 0.05, 0.10, 0.94, 1.10,
 0.11, 0.44, 0.08, 0.03, 0.93, 1.15),
 A = gl(4, 6), # 产生因子
 B = gl(6, 1, 24)  # 产生因子 注意
 )


 juice.aov<-aov(X~A+B,data=juice)
 > summary(juice.aov)
 bartlett.test(X~A, data=juice)


 rats<-data.frame(
 Time=c(0.31, 0.45, 0.46, 0.43, 0.82, 1.10, 0.88, 0.72, 0.43, 0.45,
 0.63, 0.76, 0.45, 0.71, 0.66, 0.62, 0.38, 0.29, 0.40, 0.23,
 0.92, 0.61, 0.49, 1.24, 0.44, 0.35, 0.31, 0.40, 0.56, 1.02,
 0.71, 0.38, 0.22, 0.21, 0.18, 0.23, 0.30, 0.37, 0.38, 0.29,
 0.23, 0.25, 0.24, 0.22, 0.30, 0.36, 0.31, 0.33),
 Toxicant=gl(3, 16, 48, labels = c("I", "II", "III")),
 Cure=gl(4, 4, 48, labels = c("A", "B", "C", "D"))
 )


  op<-par(mfrow=c(1, 2))
   plot(Time~Toxicant+Cure, data=rats)
    with(rats,interaction.plot(Toxicant,Cure,Time,trace.label = "Cure"))
    
    with(rats,
 interaction.plot(Cure, Toxicant, Time, trace.label="Toxicant"))
 rats.aov<-aov(Time~Toxicant*Cure, data=rats)
 bartlett.test(Time~Toxicant, data=rats)




 > ancova(Weight_Increment ~ Weight_Initial+feed ,
 data=data_feed) # 协方差分析  注意和下面的不一样，


  ancova(Weight_Increment ~ Weight_Initial*feed ,
 data=data_feed) # 协方差分析  注意和上面的不一样，


 cor.text(x, y,
 alternative=c("two.sided", "less", "greater"),
 method=c("pearson", "kendall", "spearman"),
 exact=NULL, conf.level=0.95...)


 cor.test(x,y)


 lm(formula, data, subset, weights, na.action,
 method="qr", model=TRUE, x=FALSE, y=FALSE,
 qr=TRUE, singular.OK＝TRUE, contrasts=NULL, offset, ...)
 confint(object, parm, level=0.95, ...)


 lm.reg<-lm(y~1+x)
 > summary(lm.reg)
 > confint(lm.reg, level=0.95)
 plot(lm.reg) 
 Cook 距离


 lm.pred<-predict(lm.reg,point,interval="prediction",level=0.95)
 lm.pred


 abline(lm.reg)
  res<-residuals(lm.reg)
 > plot(res)
 > text(12, res[12], labels=12, adj=(.05))


 exa<-data.frame(x1=200,x2=3000)
 lm.pred<-predict(lm.reg, exa, interval="prediction", level=0.95)


 step(object, scope, scale=0,
 direction=c("both", "backward", "forward",
 trace=1, keep=NULL, steps=1000, k=2, ...)
 lm.reg<-lm(y~x1+x2+x3+x4,data=blood)
 residuals( ), rstandard( )和rstudent( )  残差
 y.res<-residuals(lm.reg)
 y.rst<-rstandard(lm.reg)
 y.fit<-predict(lm.reg)
 plot(y.res~y.fit)


 lm.new_reg<-update(lm.reg, sqrt(.)~.) 模型按照方差稳定性修正
 coef(lm.new_reg) 提取回归系数的估计


 lm.influence(model, do.coef=TRUE) 影响函数
 ff<-lm.influence(lm.reg,do.coef=TRUE)


 cooks.distance(model, infl=im.influence(model, do.coef=FALSE),
 res=weighted.residuals(model),
 sd=sqrt(deviance(model)/df.residual(model)),
 hat=infl$hat, ...)
 dffits(model, infl=..., res=...)


 covratio(model, infl=lm.influence(model, do.coef=FALSE),
 res=weighted.residuals(model))
 influence.measures(model)


 eigen(x, symmetric, only.values=FALSE, EISPACK=FALSE)
 kappa(x, exact=FALSE, ...)
 vif(lmobj, digits=5)
 vif(lm.reg,digits=3) #共线性诊断
 cor(x2, x4)


 log<-glm(formula, family=family.generator,
 data=data.frame)


 log<-glm(formula, family = gaussian(link = identity),
 data = data.frame)


 log<-glm(formula, family = binominal(link = logit),
 data = data.frame)  #基于二项分布族的glm( )的调用格式


 accident<-data.frame(x1,x2,x3,y)
 > log.glm<-glm(y~x1+x2+x3,family=binomial,data=accident)  # logistic 回归，默认 Link=logit
 > summary(log.glm)
 log.step<-step(log.glm)
 > summary(log.step)


 log.pre<-predict(log.step, data.frame(x1=1))
 > p1<-exp(log.pre)/(1+exp(log.pre));p1


 princomp(formula, data = NULL, subset, na.action, ...) # 主成分
  student.pr<-princomp(student, cor=TRUE)
  
  factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,
 factanal( )的调用格式
 subset, na.action, start = NULL,
 scores = c("none", "regression", "Bartlett"),
 rotation = "varimax", control = NULL, ...)


 factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,
 subset, na.action, start = NULL,
 scores = c("none", "regression", "Bartlett"),
 rotation = "varimax", control = NULL, ...)


 >student<-read.table("D:/Rdata/student.txt")
 > names(student)=c("math", "phi", "chem", "lit", "his", "eng")
 > fa<-factanal(student, factors=2)


 factanal(ma,factors=3)
      # 判别分析法 马氏距离 Fisher判别法
# lda(formula, data, ... , subset, na.action) 实现Fisher 判别


iris.lda <- lda(Species ~ Sepal.Length + Sepal.Width
 + Petal.Length + Petal.Width)
 > iris.lda
 > iris.pred=predict(iris.lda) $ class


 z <- lda(group~x1+x2+x3+x4, data=w, prior=c(1, 1)/2)
  predict(z, newdata=newdata)
  dimnames(newdata)<-list(NULL, c("x1", "x2", "x3", "x4"))
  
  dist(x, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)  距离计算
  
  dim(x)<-c(5, 1)  下标转换
  > d<-dist(x)
  hc2<-hclust(d,"complete")
  
  iris.hc<-hclust(dist(iris[,1:4]))
 > plot(iris.hc,hang=1)


 plot(iris.hc,labels=FALSE,hang=-1)
 plclust(iris.hc1,labels = FALSE, hang=-1)
 re<-rect.hclust(iris.hc1,k=3)   按给定的类的个数(或阈值)进行聚类
 iris.id <- cutree(iris.hc1,3)  输出编制成若干组 
  
 cancor(x, y, xcenter = TRUE, ycenter = TRUE) # 典型相关分析


 matplot(x, z.df, ylim=c(0,4), xlab="x", ylab="",
 +         col=1:6, type="l", lwd=2)
 > text(0.20,3,"y=0")




 > 




 bioassay.post<-function(alpha=0.1,beta=5){
 +      k<-4
 +      x <- c(-0.86, -0.30, -0.05, 0.73)
 +      n <- c(5, 5, 5, 5)
 +     y <- c(0, 1, 3, 5)
 +    prod<-1
 +   prod <- prod((exp(alpha+beta*x)/(1+exp(alpha+beta*x[i])))^y
 +                    *(1/(1+exp(alpha+beta*x)))^(n-y))
 +     return(prod)}


 mlpost<-function(alpha=0.1,beta=5){-log(bioassay.post(alpha,beta))}
 > mle(mlpost)  贝叶斯估计
 for(j in 1:n){
 ra[j]<-sample(alphax,1,replace=T,prob=w)
 postb<-bioassay.post(ra[j],betay)
 wb<-postb/sum(postb)
 rb[j]<-sample(betay,1,replace=T,prob=w)
 }
 eth.tauy<-matrix(rep(0,24008),8,3001)
 > sdth.tauy<-matrix(rep(0,24008),8,3001)
 > for (j in (1:8)){
 > for (i in (1:3001)){
 eth.tauy[j,i]<-(y[j]/v[j]+muhat[i]/tausq[i])
 /(1/v[j]+1/tausq[i])
 sdth.tauy[j,i]<-sqrt(((1/tausq[i])/(1/v[j]+1/tausq[i]))
 *((1/tausq[i])/(1/v[j]+ 1/tausq[i]))
 *vmu[i]+1/(1/v[j]+1/tausq[i]))
 }
 }
 > par(mfrow=c(1,2))
 > taux<-matrix(rep(tau,8),8,byrow=T)
 > matplot(t(taux),t(eth.tauy), ylim=(c(-5,30)),
 type="l", xlab="tau",lty = 1:8, lwd = 1,col=1,
 ylab="Estimate Treatment Effects",
 main="Conditional posterior mean")
 > School<-c("A","B","C","D","E","F","G","H")
 > text(x=rep(20,8),y=t(eth.tauy)[2400,],School)
 > matplot(t(taux),t(sdth.tauy), ylim=(c(0,20)),  分层贝叶斯统计分析
 type="l", xlab="tau",lty = 1:8, lwd = 1, col=1,
 ylab="Posterior Standard Deviations",
 main="Conditional posterior SD")
 > text(x=rep(20,8),y=t(sdth.tauy)[2400,],School)