# [R experiment. 5] R program design

The solution is not single. The following method has xuanzi's personal preference, so it is only for reference. If there are any mistakes, please correct them in the comments section!

5.1 write an R program (function). Input an integer n, if n < = 0, then terminate the operation, and output a sentence: "it is required to input a positive integer"; Otherwise, if n is even, n is divided by 2 and assigned to n; Otherwise, assign 3n+1 to N, and continue to cycle until n=1, then stop the operation, and output a sentence: "the operation is successful". This example is to verify a simple theorem in number theory.

```fo <- function(n){
if(n<=0) print('A positive integer is required')
else{repeat{
if(n==1) break
else if(n%%2==0) n <- n/2
else n <- 3*n+1
}
print('The operation is successful')
}
}
> fo(0)
 "A positive integer is required"
> fo(1)
 "The operation is successful"
> fo(7)
 "The operation is successful"
```

5.2 Chevalier, a Frenchman, was once an extremely keen and powerful gambler. He especially liked two ways of gambling. In the first, roll a dice four times and bet that there will be at least one 6; For the second, roll two dice 24 times and bet that there will be at least one double 6. Note that the first bet is "good" for him: there is a more than 50% chance that 6 will appear at least once. He thinks that the second way of gambling is also beneficial to him.

(1) Write the code of a function called fourdrops(). If you roll a dice four times, at least one time is 6 points, the function will return 1, otherwise it will return 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)
0
> fourhrows(4)
1
> fourhrows(4)
 0
```

(2) Given the code of a function named twentyfourthrows(), if two dice are rolled 24 times and at least one double six appears, the function will return 1, otherwise it will return 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)
0
> twentyfourthrows(24)
0
> twentyfourthrows(24)
 1
```

(3) Write code for a function called meresix() to verify Chevalier's intuition.
Note: the first two functions should be used and a normal parameter nsim should be used to specify the number of repetitions of the gambling.

```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)
 0.25 #Many experiments have proved that the first one is more beneficial to Chevalier
> meresix2(24)
0.08333333
```

5.3 write a function called biroot() to calculate the root of a second-order polynomial, that is, find the root X of the equation ax2 +bx+c=0. Note that x may be a complex root.

```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('The solution is arbitrary')
else print('unsolvable')
}
else{
if(deta>=0) list(x1,x2)
else list(x3,x4)
}
}
```

5.4 write the program of Newton method for solving nonlinear equations, and use this program to solve nonlinear equations
For the solution of x 12 + x 22 - 5 = 0 (x 1 + 1) x 2 - (3 x 1 + 1) = 0, the initial point is (0,1), and the accuracy is required to be 10-5
Newtons<-function(fun,x,ep=1e-5,it_max = 100) ############################_ Max is the maximum number of iterations

```{
index<-0 ##Indicates whether the iteration is completed successfully and meets the accuracy requirements
k<-1 ##Iterations
while(k<=it_max)
{
x1<-x;obj<-fun(x);
x<-x-solve(obj\$J,obj\$f) ##obj\$J is the jacobj matrix of the equation (system)
norm<-sqrt((x-x1)%*%(x-x1)) ##The difference of precision between x(k+1)y and 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^2+x^2-5,(x+1)*x-(3*x+1)) ##Equations
J<-matrix(c(2*x,2*x,x-3,x+1),nrow=2,byrow=T) ##jacobj matrix, the first derivative of x1 and x2 respectively
list(f=f,J=J)
}

Newtons(funs,c(0,1))
debug(Newtons)
```

Tags: function R Language

Posted by JsusSalv on Mon, 31 May 2021 04:45:14 +0930