【R实验.5】R程序设计

刺骨的言语ヽ痛彻心扉 2022-10-16 14:27 240阅读 0赞

*解法并不单一，下列方法带有璇子个人的偏好，因此仅供参考。如有错误，欢迎在评论区斧正！*

5.1 编写一个 R 程序（函数）。输入一个整 n，如果 n<=0，则终止运算，并输出一句话：“要求输入一个正整数”；否则，如果 n 是偶数，则将 n 除 2，并赋值给 n；否则将 3n+1 赋给 n， 不断循环，直到 n=1，才停止运算，并输出一句话：“运算成功”.这个例子是为了验证数论中 的一个简单定理。

fo <- function(n){ 
      if(n<=0) print('要求输入一个正整数')
      else{ repeat{ 
        if(n==1) break
        else if(n%%2==0) n <- n/2
        else n <- 3*n+1
        }
        print('运算成功')
        }
    }
    > fo(0)
    [1] "要求输入一个正整数"
    > fo(1)
    [1] "运算成功"
    > fo(7)
    [1] "运算成功"

5.2 法国人 Chevalier 曾是一个极其敏锐和厉害的赌徒，他特别喜欢两种赌博方式。在第一 种中，将一粒骰子掷 4 次并打赌说至少会出现一次 6；对第二种，将两粒骰子掷 24 次并打赌 说至少会出现一个双 6.注意到第一种赌博对他是“有利的”：有超过 50%的可能性 6 至少会出 现一次。他认为第二种赌博方式对他也是有利的。

（1）编写一个名为 fourhrows()的函数的代码，如果将一粒骰子掷 4 次至少有一次是 6 点， 该函数将返回 1 否则返回 0.

fourhrows <- function(n){ 
      result <- sample(1:6,n,replace=T)
      for(i in result) { 
        if(i==6) return(1)
        else return(0)
      }
    }
    > fourhrows(4)
    [1]0
    > fourhrows(4)
    [1]1
    > fourhrows(4)
    [1] 0

（2）给出一个名为 twentyfourthrows()的函数的代码，如果将两粒骰子掷 24 次至少出现一 个双 6 点，该函数将返回 1 否则返回 0.

twentyfourthrows <- function(n){ 
      p <- sample(1:6,n,replace=T)
      q <- sample(1:6,n,replace=T)
      for(i in 1:n){ 
      if(p[i]==6&q[i]==6) return(1)
      else return(0)}
    }
    > twentyfourthrows(24)
    [1]0
    > twentyfourthrows(24)
    [1]0
    > twentyfourthrows(24)
    [1] 1

（3）编写一个名为 meresix()的函数的代码，来验证 Chevalier 的直觉是否正确。  
注：应当利用前面两个函数并采用一个正规的参变量 nsim 来指定该项赌博的重复次数。

meresix1 <- function(nsim){ 
      count=0
      a <- sample(1:6,nsim,replace=T,prob=rep(1/6,6))
      for(i in 1:nsim){ 
        if(a[i]==6) count <- count+1
      }
      return(count/nsim)
    }  
    
    meresix2 <- function(nsim){ 
      count=0
      a <- sample(1:6,nsim,replace=T,prob=rep(1/6,6))
      b <- sample(1:6,nsim,replace=T,prob=rep(1/6,6))
      for(i in 1:nsim){ 
        if(a[i]==6&b[i]==6) count <- count+1
      }
      return(count/nsim)
    } 
    > meresix1(4)
    [1] 0.25 #多次试验证明，第一种更有利于Chevalier
    > meresix2(24)
    [1]0.08333333

5.3 编写一个名为 biroot()的函数，用来计算一个二阶多项式的根，即求方程 ax2 +bx+c=0 的根 x.注意，x 可能是复根。

biroot <- function(a,b,c){ 
      deta <- b*b-4*a*c
      x1 <- (-b+sqrt(deta))/(2*a)
      x2 <- (-b-sqrt(deta))/(2*a)
      x3 <- (-b+sqrt(4*a*c-b*b*1i))/(2*a)
      x4 <- (-b-sqrt(4*a*c-b*b*1i))/(2*a)
      if(a==0){ 
        if(b!=0) return(-c/b)
        else if(b==0&c==0) print('解为任意值')
        else print('无解')
      }
      else{  
      if(deta>=0) list(x1,x2)
      else list(x3,x4)
      }
    }

5.4 编写求非线性方程组解的 Newton 法的程序，并用此程序求解非线性方程组  
x12 +x22 -5=0 (x1+1)x2-(3x1+1)=0 的解，取初始点为(0,1),精度要求为 10-5 .  
Newtons<-function(fun,x,ep=1e-5,it\_max=100) \#\#fun为需要求解的方程(组),x为初始解，ep为精度要求，it\_max为最大迭代次数

{ 
      index<-0 ##指示是否完成迭代成功，满足精度要求
      k<-1 ##迭代次数
      while(k<=it_max)
      { 
        x1<-x;obj<-fun(x);
        x<-x-solve(obj$J,obj$f) ##obj$J即为方程(组)的jacobj矩阵
        norm<-sqrt((x-x1)%*%(x-x1)) ##x(k+1)y与x(k)精度之差
        if(norm<ep)
        { 
          index<-1
          break
        }
        k<-k+1
      }
      obj<-fun(x);
      list(root=x,it=k,index=index,FunVal=obj$f,Jacobi=obj$J)
    }
    
    funs<-function(x)
    { 
      f<-c(x[1]^2+x[2]^2-5,(x[1]+1)*x[2]-(3*x[1]+1)) ##方程组
      J<-matrix(c(2*x[1],2*x[2],x[2]-3,x[1]+1),nrow=2,byrow=T) ##jacobj矩阵，分别对x1和x2求一阶导
      list(f=f,J=J)
    }
    
    Newtons(funs,c(0,1))
    debug(Newtons)