R语言基础

数据保存与恢复

x <- c(22,23,23,24)  y <- mean(x)#开环境：开始存储变量的数据sink('./lis/mean_test.lis')#将x，y存入mean_test.lisx   y#结束输出重定向到文件中sink()#代码续行:x为向量[11,22,3388]x <- c(11,22,33+           +22+           +3333)#向量叠加assign("x",c(11,22,55))  c(x,111,x,222) -> y

math函数

x <- c(11,2,3388)  cos(x)  sin(x)   sum(x)  mean(x)  sqrt(x)  length(x)  sort(x)

序列seq

seq(1,5)#1,2,3,4,5  seq(1,5,2)#,1,3,5  seq(from=1,to=5,by=2)#1,3,5  seq(from=1,to=5,by=2)#1,3,5  seq(from=1,to=15,length.out = 3)#1,8,15   seq(from=2,along.with = c(1:5))#2,3,4,5,6   #along.with参数中序列的长度作为要产生序列的长度，序列内容对结果无影响seq(from=2,by=2,along.with = c(1,2,5,8))#2,4,6,8    seq(from=2,by=2,along.with = c(1,2,3,8))#2,4,6,8  #序列中元素重复后拼接x<-c(9,5,6)  rep(x,2)#9,5,6,9,5,6  rep(x,times = 3)#9,5,6,9,5,6,9,5,6  rep(x,each=2)#9,9,5,5,6,6  #->在控制台不显示结果的赋值#<-在控制台显示结果的赋值c(4,5,2)->x  x>4 -> y#FALSE TRUE FALSE  #无效值或缺失值NA、NaNx<-c(1:4,NA,2:3)  is.na(x)  #字符串向量z <- c("qq","zz")#paste()接收任意参数，依次连接到字符串向量#seq代表相隔字符，默认单个空格paste(1:4)#"1" "2" "3" "4"  paste("a",1:6,seq="")#"a 1 " "a 2 " "a 3 " "a 4 " "a 5 " "a 6 "  paste("a",1:6)#"a 1" "a 2" "a 3" "a 4" "a 5" "a 6"   paste(c("a","b"),1:4,seq="")#"a 1 " "b 2 " "a 3 " "b 4 "   paste("today",date())#"today Sun Nov 20 19:25:58 2016"

索引向量

x<-c(11,22,33)  x[c(1,2,3,2,1)]#11 22 33 22 11  x[1:2]#11 22  x[-(1:2)]#33  names(x)<-c("一","二","三")#设置字符串下标索引  x[c("一")]#11  

mode为对象集

#有numeric、complex、logical、character和rawx<-c(11,22,3388)#11   22 3388  mode(x)#numeric  length(x)  as.character(x)#"11"   "22"   "3388"  a<-c(1.0-5i,29+88i)  mode(a)#complex  

设置对象属性

h=c(5:12)  attr(h,"name")<- "hello"#分配h的name属性值为hello

因子和有序因子

my_num <- c(11,22,34,71,14,58,21,22)  factor(my_num) -> nums  #[1] 11 22 34 71 14 58 21 22  #Levels: 11 14 21 22 34 58 71  levels(nums)#"11" "14" "21" "22" "34" "58" "71"  ordered(nums)  #[1] 11 22 34 71 14 58 21 22  #Levels: 11 < 14 < 21 < 22 < 34 < 58 < 71  age <- c(25,12,15,12,25)  ordered(age,levels=c(25,12))  #[1] 25   12   <NA> 12   25    #Levels: 25 < 12  score <-c(88,85,75,97,92,77,74,70,63,97)  cut(score,breaks=3)#分成三组  #[1] (85.7,97]   (74.3,85.7] (74.3,85.7]  #[4] (85.7,97]   (85.7,97]   (74.3,85.7]  #[7] (63,74.3]   (63,74.3]   (63,74.3]    #[10] (85.7,97]    #Levels: (63,74.3] (74.3,85.7] (85.7,97]  

循环语句

#forz<-c()  x<-c(1:10)  y<-c(11:20)  for (i in 1:length(x)){    z[i]=x[i]^2+y[i]^2  }  z#122 148 178 212 250 292 338 388 442 500  #whilei = 1  while(x[i]^2<10){    i=i+1    x[i]=x[i]^2  }  x#1  4  3  4  5  6  7  8  9 10

条件语句

z<-c()  x<-c(1:10)  y<-c(11:20)  for(i in 1:length(x)){    if(x[i]^3>y[i]^2){      z[i]=x[i]^3    }else{      z[i]=y[i]^2    }  }  z#121  144  169  196  225  256  343  512  729 1000

R语言科学计算

分类（组）统计

#1、准备分组数据fruit_class<-c("苹果","梨子","橘子","草莓","苹果","橘子","橘子","草莓","橘子","草莓")  fruit_prices<-c(3.5,2.5,1.5,5.5,4.2,3.2,2.8,4.8,2.9,5.8)#2、平均价格统计tapply(fruit_prices, fruit_class, mean)#3、最低价格统计tapply(fruit_prices, fruit_class, min)#4、最高价格统计tapply(fruit_prices, fruit_class, max)#5、标准差统计tapply(fruit_prices, fruit_class, sd)#标准误（描述抽样误差）：S/sqrt(n):S为样本的标准差strerr <- function(x) sqrt(var(x)/length(x))  tapply(fruit_prices, fruit_class, strerr)

数组与矩阵基础

#1、数组与矩阵的维数my_num<-c(1:9)#1 2 3 4 5 6 7 8 9  dim(my_num) <-c(3,3)  #my_num  #   ````[,1] [,2] [,3]  # [1,]    1    4    7  # [2,]    2    5    8  # [3,]    3    6    9  ###2、切片c(my_num[1,2],my_num[2,3])# 4 8  ###3、索引向量x <- array(10:20,dim=c(2,5))  #       [,1] [,2] [,3] [,4] [,5]  # [1,]   10   12   14   16   18  # [2,]   11   13   15   17   19  i <- array(c(1:3,5:4,3:5),dim=c(2,3))  #       [,1] [,2] [,3]  # [1,]    1    3    4  # [2,]    2    5    3x[i]#10 11 12 14 13 12x[i]<-111# [,1] [,2] [,3] [,4] [,5]# [1,]  111  111  111   16   18# [2,]  111  111   15   17   19#4、array函数:根据dim对向量生成数组，列优先c(1:20)->hmya<-array(h,dim=c(4,5))# 列优先 [,1] [,2] [,3] [,4] [,5]# [1,]    1    5    9   13   17# [2,]    2    6   10   14   18# [3,]    3    7   11   15   19# [4,]    4    8   12   16   20#5、数组转换为向量x<-array(c(1:10),dim=c(2,5))as.vector(x)#c(1:10)#6、matrix矩阵：数据，行数，列数，是否安行分配matrix(c(1:10),2,5,TRUE)#       [,1] [,2] [,3] [,4] [,5]# [1,]    1    2    3    4    5# [2,]    6    7    8    9   10#7、对角矩阵x<-diag(c(1:8))#对角线从1到8的8*8矩阵diag(x)#提取对角线元素1到8

数组运算

#1、四则运算:数组对应元素进行运算mya<-array(c(1:20),dim=c(2,10))myb<-array(c(2),dim=c(2,10))mya+mybmya*myb#2、向量连接x2<-c(101:105)x1<-c(1:10)cbind(x1,x2)#向量行换为列再连接#       x1  x2# [1,]  1 101# [2,]  2 102# [3,]  3 103# [4,]  4 104# [5,]  5 105# [6,]  6 101# [7,]  7 102# [8,]  8 103# [9,]  9 104# [10,] 10 105rbind(x1,x2)#列换为行再连接#     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]# x1    1    2    3    4    5    6    7    8    9    10# x2  101  102  103  104  105  101  102  103  104   105

矩阵运算

#1、矩阵连接x3<-matrix(c(1:4),2,2)x4<-matrix(c(101:104),2,2)cbind(x3,x4)#       [,1] [,2] [,3] [,4]#[1,]    1    3  101  103# [2,]    2    4  102  104 rbind(x3,x4)#       [,1] [,2]# [1,]    1    3# [2,]    2    4# [3,]  101  103# [4,]  102  104#2、矩阵转置:t(x),aperm(x,perm=c(1,2,3))array(c(1:6),dim=c(2,3))->myat(mya)#3、矩阵乘积%*%a=array(c(1:4),dim=c(2,2))b=array(c(5:8),dim=c(2,2))a%*%b#       [,1] [,2]# [1,]   23   31# [2,]   34   46#4、内积:矩阵内积即为矩阵乘积##向量内积a<-c(1:3)b<-c(4:6)crossprod(a,b)#32:1*4+2*5+3*6##向量外积a<-c(1:4)b<-array(c(1:4),dim=c(4,1))a%o%b#       [,1] [,2] [,3] [,4]# [1,]    1    2    3    4# [2,]    2    4    6    8# [3,]    3    6    9   12# [4,]    4    8   12   16#5、求线性方程组:a%*%x=b中x的值b=c(8:9)a=array(c(1:4),dim = c(2,2))solve(a,b)#-2.5  3.5#6、矩阵求逆:solve(x)array(c(1:4),dim=c(2,2))solve(array(c(1:4),dim=c(2,2)))#       [,1] [,2]# [1,]   -2  1.5# [2,]    1 -0.5###7、矩阵特征值求解:eigen(x,symmetric,only.values=FALSE)a<-array(c(1:16),dim=c(4,4))eigen(a)# $values# [1]  3.620937e+01 -2.209373e+00  1.599839e-15# [4]  7.166935e-16# $vectors# [,1]        [,2]       [,3]       [,4]# [1,] 0.4140028  0.82289268 -0.5477226  0.1125155# [2,] 0.4688206  0.42193991  0.7302967  0.2495210# [3,] 0.5236384  0.02098714  0.1825742 -0.8365883# [4,] 0.5784562 -0.37996563 -0.3651484  0.4745519#8、求解矩阵行列式det(array(c(1:4),dim=c(2,2)))#-2#9、奇异分解a<-array(c(1:16),dim=c(4,4))svd(a)

R语言计算实例

list:创建和读写数据

#1、创建学生数据集:list()mystudents<-list(name="students",class="101",stdt.ages=c(22,24,20),stdt.name=c("zhangsan","lisi","wangwu"))#2、读取列表mystudents#3、获取学生数据集的字段总数length(mystudents)#4#4、查看数据集中所有学生姓名和年龄c(mystudents)

data.frame:更优的list

#1、创建data.frame,存学生数据mysts<-data.frame(name=mystudents$stdt.name,age=mystudents$stdt.ages)#       name age# 1 zhangsan  22# 2     lisi  24# 3   wangwu  20#2、将年龄都增加一岁的操作步骤#attach:字段副本绑定在搜索路径中attach(mysts)age+1->mysts$agemysts#       name age# 1 zhangsan  23# 2     lisi  25# 3   wangwu  21# detach(mysts)

最小二乘法拟合

#lsift：#x:行对应情况，列对应变量#y:结果，可为矩阵#Wt：可选参数，权重向量#Intercept：是否使用截距项#Tolerance：公差将用于矩阵分解#Yname：用于响应变量名称y<-c(2,4,6,8)x<-c(1,2,3,4)lsfit(x,y)##y=2x+0# $coefficients# Intercept         X # 0                 2 

交叉因子频率分析

#1、划分数据分布区间:cut(x,3):x向量分成3个区间y<-c(11,22,13,14,11,22,31,31,31,14)cut(y,5)->cuty#2、使用table函数统计数据再每个区间出现的频率table(cuty)# 5       0       2       0       3 #3、使用hist函数:生成分布直方图bins<-seq(min(y),max(y),by=4)#生成序列hist(y,breaks=bins,col = "lightblue",axes=FALSE)axis(1,bins)#1代表x轴axis(2)#2代表y轴

向量模长计算

#1、模长函数定义 vector_length<-function(x1,x2,x3){  vlength<-sqrt(x1^2+x2^2+x3^2)  vlength}#2、计算向量模长vector_length(12,33,19)#39.92493#3、n维向量模长计算vector_length<-function(x){  temp<-0  for(i in 1:length(x)){    temp<-temp+x[i]^2  }  vlength<-sqrt(temp)  vlength}vector_length(c(11,22,33,44,55))#81.57818

欧氏距离:计数

mycount<-function(...){  temp<-0  for(i in c(...)){    temp=temp+1  }  temp}mycount(11,22,3,34,45)#5

统计分析基础

回归分析

#单变量线性回归x<-c(5,7,9,11,16,20)y<-c(1,2,3,4,7,9)plot(x,y)#散点图lsfit(x,y)#最小二乘法拟合# $coefficients#Intercept(截距)          X（斜率）# -1.8016529        0.5413223 # $residuals（残差）# [1]  0.09504132  0.01239669 -0.07024793 -0.15289256# [5]  0.14049587 -0.02479339abline(lsfit(x,y))#绘制散点图和回归线#lm:更详细的回归分析lm(y~x)->xysummary(xy)#Estimate Std. Error    t value   Pr(>|t|):对斜率和截距#估计值   估计标准差  假设检验t值  是否接受假设检验#多元线性回归:lm函数y=b0+b1x1+b2x2+...+bkxk+ey<-c(5,7,9,11,16,20)x<-c(1,2,3,4,7,9)x2<-c(6,8,10,13,16,20)lm(y~x+x2)->xy2summary(xy2)# Coefficients:#   Estimate Std. Error t value Pr(>|t|)# (Intercept)  2.13846    0.30735   6.958 0.006091# x            1.43077    0.10390  13.770 0.000829# x2           0.24615    0.06111   4.028 0.027500#非线性回归:非线性回归可通过变量变换变化为线性模型#y=b0+b1x+b2x^2+...+bkx^k+e#y=ae^(bx)等等非线性函数x<-c(1,2,3,4,7,8,9)y<-100+10*exp(x/2)+rnorm(x)nlmod<-nls(y~Const+A*exp(B*x))summary(nlmod)# Estimate Std. Error t value Pr(>|t|)    # Const 99.265297   0.838620  118.37 3.06e-08 ***#   A      9.996232   0.175355   57.01 5.67e-07 ***#   B      0.500196   0.001932  258.89 1.34e-09 ***# plot(x,y,main="nls(o)")#main为标题# curve(100+10*exp(x/2),col=4,add=TRUE)#画实际方程# lines(x,predict(nlmod),col=2,type = 'b')#画预测的回归线

数据分析

#1、区间频率分布mag<-c(1.6,0.9,2.1,2.2,2.3,1.7,1.3,1.6,4.7,1.2,0.9,4.7,0.6,5.3,1.1,4.8,4,4.2,4.6,1.3,2.1,1.5,3)magefactor<-factor(cut(mag,5))#分5区间table(magefactor)#区间频率统计hist(mag,breaks=5)#绘制直方图#2、数据直方图:header表示文件头做变量名，seq是分隔符earthquake<-read.table("eqweek.csv",header = TRUE,seq=",")#3、数据散点图:hist(mag)rug(mag)#4、五分位点#最大最小中位上下四分位fivenum(mag)#0.6 1.3 2.1 4.1 5.3#5、累积分布函数mag_ecdf<-ecdf(mag)#计算累积分布# Empirical CDF # Call: ecdf(mag)# x[1:18] =    0.6,    0.9,    1.1,  ...,    4.8,    5.3# plot(mag_ecdf,do.points=FALSE,verticals=TRUE)#绘制累积分布函数#6、核密度估计:densityhist(mag,prob=TRUE)lines(density(mag))

统计分析案例

绘图

#获取文件数据jiuye<-read.table("csv/work.csv",header = TRUE,sep=",")#筛选数据集中农业的数据:grepl(字符串匹配)jyhy<-jiuye$行业名称[grepl("农业",jiuye$行业名称)]jygz<-jiuye$平均劳动报酬[grepl("农业",jiuye$行业名称)]names(jygz)<-jyhy#为数据设置变量名#绘制点图dotchart(jygz)#绘制条形图barplot(jygz,horiz = TRUE)#绘制饼图pie(jygz)#茎叶图和箱型图#茎叶图:stem(变量,scale=长度,width=绘图宽度,atom=容差)cp<-read.table("csv/product.csv",header = TRUE,sep = ",")stem(cp$单位成本.元.台.,scale=2)#个位数分成两段#箱型图:boxplotboxplot(jiuye$平均教育经费)#五分位点确定箱型图

评价值

#1、平均值mean(jiuye$平均劳动报酬)#30203.17#2、加权平均值cp<-read.table("csv/product.csv",header = TRUE,sep = ",")weighted.mean(cp$单位成本.元.台.,cp$产量.台.)#322.9051#3、数据排序sort(jiuye$平均教育经费)#42  44 143 196 230 308sort(jiuye$平均教育经费,decreasing = TRUE)#308 230 196 143  44  42#4、中位数median(jiuye$平均教育经费)#169.5#极差max(jiuye$平均教育经费)-min(jiuye$平均教育经费)#266#半极差#quantile(jiuye$平均教育经费):四分位数# 0%      25%    50%    75%   100% # 42.00  68.75 169.50 221.50 308.00 #半极差即为上下四分位数之差# IQR(jiuye$平均教育经费)#152.75#方差# var(jiuye$平均教育经费)#11153.5#标准差# sd(jiuye$平均教育经费)#105.6101#变异系数：也称离散系数:标准差除以平均值

评价系数

#概率分布离散程度的归一化量度# sd(jiuye$平均教育经费)/mean(jiuye$平均教育经费)#0.6580071#样本平方和：样本校正平方和：样本与均值差的平方求和# sum((jiuye$平均教育经费-mean(jiuye$平均教育经费))^2)#55767.5#偏度系数:衡量实数随机变量概率分布的不对称性# mymean<-mean(jiuye$平均教育经费)# mysd<-sd(jiuye$平均教育经费)# myn<-length(jiuye$平均教育经费)# x<-jiuye$平均教育经费# myn/((myn-1)*(myn-2))*sum((x-mymean)^3)/mysd^3#0.08632698#峰度系数:衡量实数随机变量概率分布的峰态#正态分布时，峰度系数为0# mymean<-mean(jiuye$平均教育经费)# mysd<-sd(jiuye$平均教育经费)# myn<-length(jiuye$平均教育经费)# x<-jiuye$平均教育经费# ((myn*(myn+1))/((myn-1)*(myn-2)*(myn-3))*sum(x-mymean)^4)/mysd^4-(3*(myn-1)^2)/((myn-2)*(myn-3))#-6.25##峰度系数越接近0，数据越接近正态分布！！！！# cp<-read.table("csv/product.csv",header = TRUE,sep = ",")# mymean<-mean(cp$序号)# mysd<-sd(cp$序号)# myn<-length(cp$序号)# x<-cp$序号# ((myn*(myn+1))/((myn-1)*(myn-2)*(myn-3))*sum(x-mymean)^4)/mysd^4-(3*(myn-1)^2)/((myn-2)*(myn-3))##累积分布概率:pnorm# mymean<-mean(cp$产量.台.)# mysd<-sd(cp$产量.台.)# myn<-length(cp$产量.台.)# x<-cp$产量.台.# pnorm(x,mymean,mysd)# plot(x,pnorm(x,mymean,mysd))##概率密度函数:dnorm(变量，平均值，标准差)# mymean<-mean(cp$产量.台.)# mysd<-sd(cp$产量.台.)# myn<-length(cp$产量.台.)# x<-cp$产量.台.# dnorm(x,mymean,mysd)# plot(x,dnorm(x,mymean,mysd))##正态分布随机数:rnorm(长度，平均值，标准差)# rx<-rnorm(50,0,1)# plot(rx,dnorm(rx))##分位点:qnorm(0.25,mean,sd)# mymean<-mean(cp$产量.台.)# mysd<-sd(cp$产量.台.)# qnorm(0.25,mean=mymean,sd=mysd)#下25%分位点值# plot(x,dnorm(x,mymean,mysd))# abline(v=qnorm(0.25,mean=mymean,sd=mysd))#分位点划线##频率直方图:hist# hist(jiuye$平均劳动报酬,freq = TRUE)##经验累积分布于正态分布# plot(ecdf(jiuye$平均劳动报酬),verticals = TRUE,do.p=FALSE)# lines(x,pnorm(x,mean(jiuye$平均劳动报酬),sd(jiuye$平均劳动报酬)),col="blue")##正太检验与分布拟合#QQ图：测试数据分布是否近似某种分布，正态分布近似于直线# qqnorm(jiuye$平均劳动报酬)# qqline(jiuye$平均劳动报酬)#正态检验：W检验shapiro.test(),ks检验# shapiro.test(cp$产量.台.)#W = 0.8866, p-value = 0.3008:p值大于0.05，故认为产量正态分布# ks.test(rnorm(80),rnorm(40))#D = 0.15, p-value = 0.5725:p值大于0.05，可认为两分布为同一分布# alternative hypothesis: two-sided##其他还有很多分布在R语言中均有函数与其对应

模拟退火包使用

library("stats")#读入csv数据e<-read.csv("C:\\Users\\jerry123\\Documents\\R\\MyTest\\mle\\data.csv",header=F)e<-e[1:1655,1]T<-length(e)BIGINT<-1000000gamma<-1000G<- function(gamma,c,t) (1+exp(-gamma*(t/T-c)))^(-1)loglik<-function(param){  a0<-param[1]  deltaxing<-param[2]  c<-param[3]  ll=0  for(t in 1:T){    ll<-ll-1/2*log(a0+deltaxing*G(gamma,c,t))-1/2*e[t]^2/(a0+deltaxing*G(gamma,c,t))    ll<-ll+(min(a0+deltaxing*G(gamma,c,t),0))^2  }  ll<-ll+BIGINT*((min(1-c,0))^2+(min(c,0))^2)  return(ll)}set.seed(123)res<-optim(par=c(2,2,0.8),fn=loglik,method = "SANN",control = list(maxit = 200))summary(res)