R语言函数编写

语言函数编写是一个非常重要的技能，它能在没有任何包能实现自己的想要的功能的情况下，实现问题的解决。
在R中，不同行的语句会被视为不同的命令，如果想在同一行内放置多个命令函数，则需要在每个函数之间加上 ; 进行分割

控制流

R函数中的控制流函数和perl语言等都十分相像，在格式上有一些差别。

• if
  if函数：if(条件){语句}else if(条件){语句}else{语句}
  在if函数中，只要else if中有一个条件满足后，后面所有的else if都会直接掠过。
• repeat函数
  repeat不能自己跳出循环，要在使用了break命令后才能实现循环的结束

a=1
repeat{
  a=a+1
  if(a>10){
    print(a)
    break
}}

使用循环函数编写一个Fibonacci（斐波那契数列）数的程序：

a<-c(1,1);
repeat{
  if(a[length(a)]>1000){
    a<-a[-length(a)];break
  }
  a[length(a)+1]<-a[length(a)-1]+a[length(a)];
  }

>a
[1]   1   1   2   3   5   8  13  21  34  55  89 144 233 377 610 987

自定义函数

在R中使用function函数进行自定义函数的编写。

• Usage
  function_name<-function( arglist ) expr
  return(value)
• Arguments
  arglist:Empty or one or more name or name=expression terms.
  expr:An expression.
  value:An expression.

定义一个计算x*y+z结果的函数：

> fun<-function(x,y,z){
  a=x*y+z
  return(a) #在R中，可以不写return函数，在没有return函数的时候，会默认返回最后函数最后一行的值
}
> fun(1,2,3)
[1] 5

编写一个使用二分法求非线性方程根的函数

image.png
二分法：
对于区间[a，b]上连续不断且f（a）·f（b）<0的函数y=f（x），通过不断地把函数f（x）的零点所在的区间一分为二，使区间的两个端点逐步逼近零点，进而得到零点近似值的方法叫二分法。
二分法求解非线性函数在区间[1,2]上的根的R脚本：

fun<-function(x){
  x^3-x-1;
}
dichotomy<-function(fun,a,b,eps=1e-6){
  if(fun(a)*fun(b)>0){ 
    list(fail="寻找区间内没有根")
  }else{
    while((b-a)>1e-6){#精度要达到1e-6
      if(fun((a+b)/2)*fun(a)<0){
        b=(a+b)/2;
      }else{
        a=(a+b)/2
      }
    }
  }
  return((a+b)/2)
}

> dichotomy(fun,1,2)
[1] 1.324718

使用牛顿迭代公式求解非线性方程组

牛顿迭代算法原理： 林业科学研究院-符利勇老师-R语言与林业统计建模课程课件

https://www.matongxue.com/madocs/205.html # 马同学高等数学 如何通俗易懂地讲解牛顿迭代法？
https://www.matongxue.com/madocs/7.html # 马同学高等数学 如何通俗地解释泰勒公式？
https://www.guokr.com/question/461510/ #牛顿开方法的算法及其原理
https://baike.baidu.com/item/%E9%9B%85%E5%8F%AF%E6%AF%94%E7%9F%A9%E9%98%B5/10753754 #jacobian矩阵
http://jacoxu.com/jacobian%E7%9F%A9%E9%98%B5%E5%92%8Chessian%E7%9F%A9%E9%98%B5/ #Jacobian矩阵和Hessian矩阵
https://wenku.baidu.com/view/cac8012db4daa58da0114ae8.html #二元非线性方程组求根的牛顿迭代法

牛顿迭代法和二分法有一定的相似性，它们都是通过不断的迭代使得得到的值，不断趋近于真实阶。在牛顿迭代法中，是不断的得到函数图形的切线而来缩短与真实解的距离。
实现—课件内容

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,ncol=2,byrow=T)
  return(list(f=f, J=J))
}
Newtons<-function(funs,x,eps=1e-5,it_max=100){
  index<-0
  k<-1
  while(k<=it_max){
    x1<-x
    f.value<-funs(x)
    x<-x-solve(f.value$J,f.value$f)
    f.e<-sqrt((x-x1)%*%(x-x1))
    if (f.e<eps){
      index<-1
      break}
    k<-k+1}
  f.value<-funs(x)
  return(list(root=x,it_max=k,index=index,f.value=f.value$f))}

> Newtons (funs,c(0,1))
$`root`
[1] 1 2

$it_max
[1] 6

$index
[1] 1

$f.value
[1] 1.598721e-14 6.217249e-15