R实验.5R程序设计

Posted 2021-06-07 是璇子鸭

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

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)