20160921第二次非参作业

P318 A.5 利用二分法来求解方程

f <- function(x){
  #待求解的方程
  f_x<- 2*x^3-4*x^2+3*x-6 
  return(f_x)
}
solve <- function(a,b,eps=1e-06){
  #x,y为输入的求解区间,默认精度为1e-06
  if(f(a)*f(b)>0) return("error")
  while(abs(a-b)>=eps)
  {
    x=(a+b)/2
    if(f(x)==0){break}
    else if(f(a)*f(x)<0) b=x
    else a=x
  }
  return(x)
}
solve(-10,10)#求解方程

## [1] 2

P318 A.7 编码问题

###编码函数
coding <- function(string_sample)
{
  library(stringr)
  library(stringi)
  length_string <- str_length(string_sample) #明文字符串长度
  a=NULL#初始化储存字符的对象
  for(i in 1:length_string)#将字符串分割成一个个字符
  {
    a[i]=stri_sub(string_sample,i,length=1)
  }
  for(i in 1:length_string)#按照规则进行编码
    {
      for(j in 1:26)
        {
            if(j>=1 && j<=13){
                if(a[i]==letters[j]) {
                      a[i]=LETTERS[j+13]
                      break}
                if(a[i]==LETTERS[j]){
                      a[i]=letters[j+13]
                      break}
            }else{
                if(a[i]==letters[j]){
                      a[i]=LETTERS[j-13]
                      break}
                if(a[i]==LETTERS[j]){
                      a[i]=letters[j-13]
                      break}
              }
        }
    }
  string_result=NULL #初始化储存最终结果的对象
  for(i in 1:length_string)
    {
        string_result<- paste(string_result,a[i],sep = '')
    }
  print(string_result)
}
###演示编码
coding("People")

## [1] "cRBCYR"

P319 A.9 报号退圈

b=rep(1,13)#这13个人是否在圈里，在圈里为1，不在圈里为0
i=0 #初始化“报号次数总计”
while(sum(b)>1)#如果圈子里的人多于1人的话，继续报号
{
  for(j in 1:13)#十三个“位置”一循环
  {
    if(b[j]!=0){#如果位置上没有人(即b[j]=0)的话则跳过
    #如果位置上有人(即b[j]=1)的话继续
    i=i+1 #报号
    if(i%%3==0) b[j]=0 #报到3者退出圈子
    }
  }
}
a=1:13 #13个人的相对位置
a[b!=0]#最终留在圈里的人

## [1] 13

P319 A.11

###必要准备

#导入数据
data_0921 <- read.table("E:\\文档\\102-大学\\003-大三\\非参数统计\\2016.09.21-第一次个人作业\\basket.txt",header = T)
#判定元素是否为A,B
is_A<-data_0921[2:6]==LETTERS[1]
is_B<-data_0921[2:6]==LETTERS[2]
#计算每行A和B的个数
count_A <- apply(is_A,1,sum)
count_B <- apply(is_B,1,sum)
#将上述计算结果合并到原始导入的数据集中
data_0921<-within(data_0921,{
  A<-count_A
  B<-count_B})

###(1).一行中至少同时出现一个A和一个B的行数
nrow(subset(data_0921,A>=1&B>=1))

## [1] 12

###(1).一行中至少同时出现一个A和一、三个B的行数
nrow(subset(data_0921,A>=1&B>=3))

## [1] 3

补充练习1 正态分布的条件分布

###数据导入和处理准备
data_0921 <- read.csv("E:\\文档\\102-大学\\003-大三\\非参数统计\\2016.09.21-第一次个人作业\\protein（补充作业1数据）.csv",header=T)
data_0921_1 <- data_0921[,2:4]
data_0921_1 <<- as.matrix(data_0921_1)

###编写四个估计的R计算函数
E_v1_v2v3 <- function(data,x2x3)#求条件期望E(V1|(V2,V3))的mle
{
  #求均值向量miu的极大似然估计
  miu_mle <- apply(data,2,mean) 
  #求协方差阵var的极大似然估计
  fun1 <- function(x)
  {
    y=(x-miu_mle)%*%t(x-miu_mle)
    return(y)
  }
  s=matrix(rep(0,9),ncol = 3)
  for(i in 1:nrow(data))
  {
    s=s+fun1(data[i,])
  }
  var_mle <- s/nrow(data)
  #求条件期望E(V1|(V2,V3))的mle
  mle_e <- miu_mle[1]+var_mle[1,2:3]%*%solve(var_mle[2:3,2:3])%*%(x2x3-miu_mle[2:3])
  return(mle_e)
}
Var_v1_v2v3 <- function(data,x2x3)#求条件方差Var(V1|(V2,V3))的mle
{
  #求均值向量miu的极大似然估计
  miu_mle <- apply(data,2,mean) 
  #求协方差阵var的极大似然估计
  fun1 <- function(x)
  {
    y=(x-miu_mle)%*%t(x-miu_mle)
    return(y)
  }
  s=matrix(rep(0,9),ncol = 3)
  for(i in 1:nrow(data))
  {
    s=s+fun1(data[i,])
  }
  var_mle <- s/nrow(data)
  #求条件方差Var(V1|(V2,V3))的mle
  mle_v <- var_mle[1,1]-var_mle[1,2:3]%*%solve(var_mle[2:3,2:3])%*%var_mle[2:3,1]
  return(mle_v)
}
E_v2v3_v1 <- function(data,x1)#求条件期望E((V2,V3)|V1)的mle
{
  #求均值向量miu的极大似然估计
  miu_mle <- apply(data,2,mean) 
  #求协方差阵var的极大似然估计
  fun1 <- function(x)
  {
    y=(x-miu_mle)%*%t(x-miu_mle)
    return(y)
  }
  s=matrix(rep(0,9),ncol = 3)
  for(i in 1:nrow(data))
  {
    s=s+fun1(data[i,])
  }
  var_mle <- s/nrow(data)
  #求条件期望E((V2,V3)|V1)的mle
  mle_e <- miu_mle[2:3]+var_mle[2:3,1]%*%solve(var_mle[1,1])%*%(x1-miu_mle[1])
  return(mle_e)
}
Var_v2v3_v1 <- function(data,x1)#求条件方差Var((v2,v3)|v1)的mle
{
  #求均值向量miu的极大似然估计
  miu_mle <- apply(data,2,mean) 
  #求协方差阵var的极大似然估计
  fun1 <- function(x)
  {
    y=(x-miu_mle)%*%t(x-miu_mle)
    return(y)
  }
  s=matrix(rep(0,9),ncol = 3)
  for(i in 1:nrow(data))
  {
    s=s+fun1(data[i,])
  }
  var_mle <- s/nrow(data)
  ##求条件方差Var((v2,v3)|v1)的mle
  mle_v <- var_mle[2:3,2:3]-var_mle[2:3,1]%*%solve(var_mle[1,1])%*%var_mle[1,2:3]
  return(mle_v)
}

###使用函数
miumle_1 <- E_v1_v2v3(data=data_0921_1,x2x3=c(2,3))#求条件期望E(V1|(V2,V3))的mle
varmle_1 <- Var_v1_v2v3(data=data_0921_1,x2x3=c(2,3))#求条件方差Var(V1|(V2,V3))的mle
miumle_2 <- E_v2v3_v1(data=data_0921_1,x1=4)#求条件期望E((V2,V3)|V1)的mle
varmle_2 <- Var_v2v3_v1(data=data_0921_1,x1=4)#求条件方差Var((v2,v3)|v1)的mle

补充练习2 Euclid算法

gcd=function(a,b)
{
  #运用函数的递归调用
  maxab=max(a,b)
  minab=min(a,b)
  if(maxab%%minab!=0){
  x=maxab%%minab
  z=gcd(x,minab)}
  else z=minab
  return(z)
}
gcd(978549243672,849866359923)

## [1] 92463

20160921第二次非参作业

黄涛

2016年9月21日