函数基本格式：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))
