函数基本格式：name <- function(arg_1, arg_2, ...) expression

2.9.1 编写自己的函数

例2.4 编写一个二分求非线性方程根的函数，并求方程x3 ? x ?1 = 0在区间[1,2]内的根，精度为10^-6

1.建立函数

分析：输入参数有4个，f-函数，区间端点a,b，以及精度(当区间长度小于指定要求时，停止计算)

默认精度为10^-5

fzero<-function(f,a,b,eps=1e-5){
if(f(a)f(b)>0)
list(fail="finding root is fail!")
else{
repeat{
if(abs(b-a)<eps) break
x<-(a+b)/2
if(f(a)f(x)<0) b<-x else a<-x
}
}
list(root=(a+b)/2,fun=f(x))
}
f<-function(x) x^3-x-1
fzero(f,1,2,1e-6)

那么我们也可以偷懒下，因为R语言也提供了求一元方程根的函数uniroot()

uniroot(f, interval, ...,lower = min(interval), upper = max(interval),f.lower = f(lower, ...), f.upper = f(upper, ...),

extendInt = c("no", "yes", "downX", "upX"), check.conv = FALSE,

tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0)

uniroot(f,c(1,2))

例2：已知两个样本A,B。计算两个样本的T统计量

A: 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97 80.05 80.03 80.02 80.00 80.02

B: 80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97

统计量公式：当两个样本的方差相同且未知，

twosam <- function(y1, y2) {
n1 <- length(y1); n2 <- length(y2)
yb1 <- mean(y1); yb2 <- mean(y2)
s1 <- var(y1); s2 <- var(y2)
s <- ((n1-1)s1 + (n2-1)s2)/(n1+n2-2)
tst <- (yb1 - yb2)/sqrt(s*(1/n1 + 1/n2))
tst
}

A <- c(79.98, 80.04, 80.02, 80.04, 80.03, 80.03, 80.04,
79.97, 80.05, 80.03, 80.02, 80.00, 80.02)
B <- c(80.02, 79.94, 79.98, 79.97, 79.97, 80.03, 79.95,
79.97)

twosam(A,B)

后续会用T统计量估计是否相同

2.9.2 定义新的二元运算

二元函数形式：%anything% 设想，y是两个向量，

定义<x,y>=exp(-||x-y||^2/2) 其运算符为%！%

"%!%"<-function(x,y){exp(-0.5(x-y)%%(x-y))}
"%!%"(1,1)
"%!%"(2,1)

2.9.3 有名参数与缺省

下面是伪代码

fun1 <- function(data, data.frame, graph, limit) {
[function body omitted]
}

则下面三种调用方法等价的

ans <- fun1(d, df, TRUE, 20)
ans <- fun1(d, df, graph=TRUE, limit=20)
ans <- fun1(data=d, limit=20, graph=TRUE, data.frame=df)

例2.6 （图）

ep是精度，缺省时1e-5,it_max最大迭代次数，缺省时为100

Newtons<-function(fun,x,ep=1e-5,it_max=100){
index<-0;k<-1
while(k<=it_max){
x1<-x;obj<-fun(x);
x<-x-solve(obj$J,obj$f);
norm<-sqrt((x-x1)%*%(x-x1))
if(norm<ep){
index<-1;break
}
k<-k+1
}
obj<-fun(x);
list(root=x,it=k,index=index,FunVal=obj$f)
}

root-方程解的近似值， it-迭代次数，indexl指标，1-成功 0-失败 FunVal是方程root处的函数值。

funs<-function(x){
f<-c(x[1]^2+x[2]2-5,(x[1]+1)x[2]-(3x[1]+1))
J<-matrix(c(2x[1],2x[2],x[2]-3,x[1]-1),nrow=2,byrow=T)
list(f=f,J=J)
}

Newtons(funs,c(0,1))