统计模拟及其R实现-肖枝洪、朱强

第二章代码2.1-2.7

3**2+sqrt(4)
x1<-0:10;x1#如果一条语句后继续写语句，需要用；隔开
#R区分大小写
x1<-0:100;x1
x2<-x1*2*pi/100;x2
y<-cos(x2);y
plot(x1,y,type='l')
vect<-c(10,6,4,7,8)
mean(vect)
sd(vect)
min(vect)
median(vect)
max(vect)
boxplot(vect)#绘制数据盒型图
A<-matrix(c(1,2,7,3),ncol=2,byrow=T);A#ncol=2表示创建的矩阵A应该有2列,如果ncol参数被省略，那么R会自动计算矩阵的列数，使得矩阵元素能够按照给定的顺序填满矩阵
#byrow=T表示按行排列，如果是byrow=F则表示按列排列
Ai<-solve(A);Ai#求矩阵A的逆
b<-c(2,3);b
x=Ai%*%b;x#%*%表示矩阵的乘法
source("D:\\R\\统计模拟及R实现\\2.1.txt")
x<-c(1:3,10:13);x
x1<-c(1,2);x1
x2<-c(3,4);x2
x<-c(x1,x2);x
#x<-c(x1,x2);x
#[1] 1 2 3 4
#[1]表示改行第一个数的下标
x<-c(12:78);x
x1=length(x);x1
x<-c(1,2,6.25);x
y=x**2+1;y
(c(1,2,3)+c(10,20,30))/c(2,4,5)#等长度向量间进行加减乘除和乘方运算，其含义是对应元素间进行四则运算
5%/%3#%/%表示整除
5%%3#%%表示求余数
c(1,2,3)+c(10,20,30,-2,-4,-9,0,0,1)#长度不同，但长向量长度为短向量整数倍，短向量循环使用
x<-c(1,2,6.25)
sqrt(x)#sqrt、log、exp、sin、cos、tan都可以直接对向量进行运算
#options(digits=3)可以控制结果的有效位数输出
options(digits=3)
sqrt(x)
#min和max分别返回向量最小值和最大值，sum计算元素和，mean计算平均值，var计算样本方差，sd计算标准差
#若x为矩阵，var(x)为协方差阵，range返回包含最大值和最小值的向量，sort(x)返回按x的元素从小到大排序的结果向量
#order(x)返回使得x元素从小到大排序的元素的下标向量，x[order(x)]等效于sort(x);
z<-c(-3,-9,8,2,0,-8);z
range(z)
sort(z)
order(z)
z[order(z)]
sort(z,decreasing=TRUE)#从大到小排序
c(1,2,NA)-c(1,NA,4)#任何函数值与缺失值的运算结果仍为缺失值
numeric(9)#numeric(n)产生填充了n个0的向量
n=5
1:n
n:2
1:n-1
1:(n-1)
seq(-2,3)
seq(from=-2,to=3)
seq(0,1,0.1)
seq(0,by=0.2,1)
seq(by=0.2,from=0,1)
seq(0,length=5)
seq(length=5,to=8)
w <- c(4,-2,7,0)
seq(along=w)
rep(w,2)
rep(c(10,28),c(2,3))
l <- c(TRUE,TRUE,FALSE);l
l1<-w>3;l1
log(10*w)
log(10*w)>w
(w>3)&(c(0,3,5,6)<2)
!(w>3)
all(!(w>3))#判断一个逻辑向量是否都为真值
any(!(w>3))#判断一个逻辑向量是否有真值
is.na(c(1,0,NA))#is.na(x)判断x的函数是否有缺失
age=30
c("young","old")[(age>65)+1]
c1=c("x","tan(x)","年龄");c1#向量元素可以去字符串值
paste("I","love","homeland")
paste(c("X","Y"),"=",1:4)
paste("result.",1:5,seq="")
paste(c("123","234"),collapse="*")
x=seq(1,8,by=2);x[2]
(c(2,4,6)+(-2))[2]
x[2]=125;x
x[c(1,3)]
x[1:2]
x[c(1,3,2,1)]
c("a","b","c")[rep(c(2,1,3),3)]
x[c(-1,-3)]
x>3
x[x>3]
x[x<(-9)]
ages=c(Li=33,Zhang=29,Liu=18);ages
ages["Zhang"]
ages[c("Zhang","Li")]
x[c(1,3)]=c(188,189)
x
x[c(1,3)]<-0;x
x[]=0;x
x<-0;x
y=numeric(length(x))
y[x<0]=1-x[x<0]
y[x>=0]<-1+x[x>=0]
x<-c(-2.3,4,-5,7);y
matrix(data=NA,nrow=1,ncol=1,byrow=FALSE,dimnames=NULL)
#ncol为列数；byrow表示按行输入还是按列输入，默认为按列
A<-matrix(1:12,ncol=4,nrow=3,byrow=TRUE);A
B<-matrix(c(1,2),ncol=3,nrow=4);B
c(A)
x1=rbind(c(1,2),c(3,4));x1
x2=cbind(c(1,2),c(3,4));x2
A[1:2,]
A[1:2,c(1,2,4)]
rownames(A)<-c("a1","a2","a3")
colnames(A)<-paste("x",1:4,seq="")
A
A["a2","x4"]#不能运行
A+c(100,200,300)#结果是矩阵A按列拉直后的向量与向量c(100,200,300)相加
t(A)
B=matrix(c(1,0),nrow=4,ncol=2,byrow=T)
A%*%B
C=matrix(1:9,nrow=3,ncol=3,byrow=T)
diag(C)
diag(c(1:3))
diag(3)
#apply(x,margin,fun,...)apply可以对每行或每列进行计算
apply(A,1,sum)
z=1:24
dim(z)<-c(2,3,4)
z
dim(z)<-24
z
array(x,dim=length(x),dimnames=NULL)
z<-array(1:24,dim=c(2,3,4))
z
#outer(a,b,f)是一个一般性的外积函数
x=seq(-2,2,length=20)
y=seq(-pi,pi,length=20)
f<-function(x,y)  cos(y)/(1+x^2)
z<-outer(x,y,f)
persp(x,y,z)
d1<-outer(0:9,0:9);d1
d2<-outer(d1,d1,"-")
fr<-table(d2)
plot(as.numeric(names(fr)),fr,type="h",xlab="Determinant",ylab="frequency")
x=c("男","女","男","男","女")
y=factor(x);y
#factor(x,levels=sort(unique(x),na.last=TRUE),labels,exclude=NA,ordered=FALSE)
f<-factor(c(1,0,1,1,0),levels=c(1,0),labels=c("男","女"))
sex=factor(c("男","女","男","男","女"))
res.tab<-table(sex);res.tab
#tapply(x,INDEX,FUN=NULL,...,simplify=TRUE)
h<-c(165,170,168,172,159)
tapply(h,sex,mean)
#gl(n,k,length=n*k,labels=1:n,orderd=FALSE)
gl(2,3)
gl(2,1,6)
rec=list(name="黎明",age=30,score=c(85,76,90));rec
rec[[2]]
rec[[3]][2]
rec[2]
rec[["age"]]
rec$name
list.A<-list(name="wangming",age=31,weight=45)
list.B<-list(school="Qianjin xiao xue",address="Wuhan",grade="class 3")
list.AB<-c(list.A,list.B)
ev=eigen((1:3)%o%(1:3));ev
y=matrix(c(1,2,3,4,5,6,7,8,9),nrow=3,ncol=3)
eigen(y)
det(y)
d=data.frame(name=c("黎明","张聪","王建"),age=c(30,35,28),height=c(180,172,177));d
cat(c("AB","c"),c("E","F"),"\n",seq="")
#cat()可以指定一个参数"file=文件名"，把结果写到指定文件夹中，若指定文件夹存在，则覆盖原来的内容，加上append=TRUE可以在原文件后边加上
#sink指定一个文件用来把后续输出转向到该文件
#scan从指定文件读入
#read.table也可以读入
{
  x<-15
  x
}
w<-matrix(-10:4,ncol=5);w
y<-numeric(length(w));y
y[w>0]<-1
y[w<=0]<-0
s<-0
for(i in seq(along=w)){
  cat('w(',i,')=',w[i],'\n',seq='')
  s<-s+w[i]
}
t<-numeric(365)
n=23
for(i in 1:365){
  t[n]<-1
  for(j in 0:(n-1)){
    t[n]<-t[n]*(365-j)/365
  }
  t[n]<-1-t[n]
}
t
system.time(process)
t<-numeric(365)
for(i in 1:365){
  t[n]<-1-prod((365:(365-n+1))/365)
}
t
x<-1-cumprod((365:1)/365)
x

2.8

#ls函数查看当前工作空间保存的变量和函数
ls()
#rm函数可以剔除不想要的对象
hello<-function(){
  cat("Hello,world! \n")
  cat("\n")
}
hello
x=list(1,"abc")
x
f=function(x){x[[2]]="!!";x}
f(x)
x<-2
f<-function(){
  print(x)
  x<-20
}
f()
x
f01<-function(x,y){
  cat('x=',x,'\n')
  cat('y=',y,'\n')
  z=x-y
  cat('z=',z,'\n')
  z
}
f01(1,2)
f01<-function(x,y){
  browser()
  z=x-y
  z
}
f01(11:12,1:3)
debug(f01)
undebug(f01)
#例2.1
open.account=function(total){
  list(
    deposit=function(amount,total){
      if(amount==0)
        stop("Deposits must be positive! \n")
      total<-total+amount
      cat(amount,"deposited. your balance is",total,"\n\n")
    },
    withdraw=function(amount,total){
      if(amount>total)
        stop("You don't have that much money! \n")
      total<-total-amount
      cat(amount,"withdraw. Your balance is",total,"\n\n")
    },
    balance=function(total){total=total
      cat("Your balance is",total,"\n\n")
      total}
  )
}
#######
total=100
ross=open.account(total);ross
robert=open.account(200);robert
ross $ withdraw (30,100)
#例2.2
#双三次核函数
kernel.dcube=function(u){
  y=numeric(length(u))
  y[abs(u)<1]=(1-abs(u[abs(u)<1])**3)**3
  y[abs(u)>=1]=0
  y
}
kernel.smooth1=function(X,Y,kernel=kernel.dcube,h=1,m=100,plot.it=T){
  x=seq(min(X),max(X),length=m)
  fx<-numeric(m)
  for(j in 1:m){
    fx[j]<-sum(kernel((x[j]-X)/h)*Y)/sum(kernel((x[j]-X)/h))
  }
  if(plot.it){
    plot(X,Y,type="p")
    lines(x,fx)
  }
  cbind(x=x,fx=fx)
}
##
x=seq(-3,3,length=100)
y=x
f=function(x,y,ssq1=1,ssq2=2,rho=0.5){
  det1=ssq1*ssq2*(1-rho^2)
  s1=sqrt(ssq1)
  s2=sqrt(ssq2)
  1/(2*pi*sqrt(det1))*exp(-0.5/det1*(ssq2*x^2+ssq1*y^2-2*rho*s1*s2*y))
}
z=outer(x,y,f)
persp(x,y,z)
contour(x,y,z)
image(x,y,z)