R语言--(7)--自定义函数

变量名 <- function(arg_1, arg_2, ...) expression

例1 求一元方程的根

编写一个用二分法求非线性方程根的函数，并求方程 $$x^3-x-1=0$$ 在区间 $[a,b]$ 内的根，精度要求 $ϵ=10^{-6}$
绘制函数图像

x=seq(-5,5,0.1);x
f=x^3-x-1
plot(x,f,type="l", lwd=2,col="red",ylim=c(-2,2),main="f=x^3-x-1")

解题思路
取初始化区间 $[a,b]$, 当 $f(a)$ 与 $f(x)$异号，作二分法计算；否则停止计算(输出计算失败信息)

二分法计算过程
取中点 $x=\frac{a+b}{2}$, 则 $f(a)$ 与 $f(x)$ 异号，则 b=x; 否则 a = x. 当区间长度小于指定要求时，停止计算.

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      # abs() 取绝对值
			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)

$root 相当于 x 的值 $fun 相当于 y 的值

R语言自带求一元方程根的函数 uniroot()

uniroot(f, interval, lower=min(interval),upper=max(interval), tol = .Machine$double.eps^0.25, maxiter = 1000, ...)

参数 f 是定义的函数

例2：计算两样本T统计量

两个分布

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)

当两样本的方差相同，且未知，则T统计量的计算公式为 $$T=\frac{(\bar{X}-\bar{Y})}{S\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}$$ 其中 $$S^2=\frac{(n_1-1)S_1^2+(n_2-1)S_2^2}{n_1+n_2-2}$$
$\bar{X},\bar{Y}$ 分别是两组数据的样本均值
$S_1^2,S_2^2$ 分别是两组数据的样本方差
$n_1,n_2$ 分别为两组数据的个数

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)
	(yb1 - yb2)/sqrt(s*(1/n1 + 1/n2))
}

调用

twosam(A,B)

自定义二元运算 %anything%

设x,y 是两个向量，定义x与y的内积
$$⟨x,y⟩=\exp (-\Vert x-y{\Vert }^2/2)$$

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

有名参数与省略

“name=object" 给出被调用函数中的参数，则这些参数可以按照任意顺序给出

fun1 <- function(data, data.frame, graph, limit){
	...boby 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)

使用 Newton 求解非线性方程组解

Newton 牛顿迭代法
方程 $$\{\begin{array}{c}{x}_1^2+x_2^2-5=0\\ (x_1+1)x_2-(3x_1+1)=0\end{array}$$的解，取初始点 $x^{(0)}=(0,1{)}^T$, 精度要求 $ϵ=10^{-5}$

解：求解非线性方程组 $$f(x)=0,f:R^n\to R^n\in C^1$$ 的Newton法的迭代格式为
$$x^{(k+1)}=x^{(k)}-[J(x^{(k)}){]}^{-1}f(x^{(k)}),k=0,1,\dots$$
其中 $J(x)$ 为函数 $f(x)$ 的 Jacobi 矩阵：
$$J(X)=\left[\begin{array}{cccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}& \dots & \frac{\partial f_1}{\partial x_n}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial 2_1}{\partial x_2}& \dots & \frac{\partial f_n}{\partial x_n}\\ \frac{\partial f_n}{\partial x_1}& \frac{\partial f_n}{\partial x_2}& \dots & \frac{\partial f_n}{\partial x_n}\end{array}\right]$$

Newtons <- function(fun, x, ep=1e-5, it_max=1000){
	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)
}

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, byrow=T);
	list(f=f,J=J);
}

Newtons(funs, c(0,1))

递归函数计算积分

使用递归函数计算数值积分 ${\int }_1^5\frac{dx}{x}$，精度要求 $ϵ=10^{-6}$
解题思路
采用自动选择步长的复化梯形公式，方法为：每次将区间二等分，在子区间上采用梯形求积公式，如果计算满足精度要求或达到最大迭代次数，则停止计算；否则继续将区间对分

area <- function(f, a, b, eps=1.0e-6, lim = 10){
	fun1 <- function(f,a,b,fa,fb, a0, eps, lim, fun){
		d <- (a+b)/2;
		h <- (b-a)/4;
		fd <- f(d)
		a1 <- h * (fa + fd);
		a2 <- h * (fd + fb)
		if (abs(a0 - a1 -a2) < eps || lim == 0)
			return(a1 + a2)
		else{
			return(fun(f,a,d,fa,fb,a1,eps,lim-1,fun) + fun(f,d,b,fd,fb,a2,eps,lim-1,fun))
		}
	}
	fa <- f(a);
	fb <- f(b);
	a0 <- ((fa + fb) * (b-a))/2
	fun1(f,a,b,fa,fb,a0,eps,lim,fun1)
}

f <- function(x){
	1/x
}

quad <- area(f,1,5); quad