函數(shù)基本格式：name <- function(arg_1, arg_2, ...) expression
2.9.1 編寫自己的函數(shù)
例2.4 編寫一個二分求非線性方程根的函數(shù)，并求方程x3 ? x ?1 = 0在區(qū)間[1,2]內(nèi)的根，精度為10^-6
1.建立函數(shù)
分析：輸入?yún)?shù)有4個，f-函數(shù)，區(qū)間端點a,b，以及精度(當區(qū)間長度小于指定要求時，停止計算)
默認精度為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語言也提供了求一元方程根的函數(shù)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統(tǒng)計量
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
統(tǒng)計量公式：當兩個樣本的方差相同且未知，
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)
后續(xù)會用T統(tǒng)計量估計是否相同
2.9.2 定義新的二元運算
二元函數(shù)形式：%anything% 設(shè)想，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 有名參數(shù)與缺省
下面是偽代碼
fun1 <- function(data, data.frame, graph, limit) {
[function body omitted]
}
則下面三種調(diào)用方法等價的
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最大迭代次數(shù)，缺省時為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-迭代次數(shù)，indexl指標，1-成功 0-失敗 FunVal是方程root處的函數(shù)值。
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))
