刚看完《Introduction to Scientific Programming and Simulation Using R》这本书，
第一部分是R编程入门； 第二部分是数值计算，主要是解方程，求积分和优化；第三部分是概率和统计，主要讲概率、随机变量等概念和参数估计； 第四部分是simulation，主要讲Monte Carlo积分和方差降低。
Monte Carlo方法由冯·诺伊曼于二战时提出，1777年法国人布丰用此思路估计pi值被认为是Monte Carlo的起源，这个方法简单又快速，通过产生随机数，将数值计算问题变成随机变量的某些特征值（比如期望值）。
积分运算，和估计pi值一样，用hit and miss法估计。
hit_miss <- function(fun, lower, upper, n=500) {
# Monte Carlo integration using the hit and miss method
x <- runif(n, lower, upper)
f.value <- sapply(seq(lower, upper, 0.001), fun)
f.min <- min(f.value)
f.max <- max(f.value)
y <- runif(n, f.min, f.max)
hits <- sum(y <= fun(x))
area <- (upper-lower) * f.min + hits/n * (upper-lower) * (f.max-f.min)
return(area)
}
hit and miss方法收敛太慢，效率并不高，通常所说的MC积分是指下面这个方法。
Optimization means to seek minima or maxima of a funtion within a given defined domain.
If a function reach its maxima or minima, the derivative at that point is approaching to 0. If we apply Newton-Raphson method for root finding to f', we can get the optimized f.
The foundamental idea of numerical integration is to estimate the area of the region in the xy-plane bounded by the graph of function f(x). The integral was esimated by dividing x into small intervals, then adds all the small approximations to give a total approximation.
Trapezoidal rule Numerical integration can be done by trapezoidal rule, simpson’s rule and quadrature rules. R has a built-in function, integrate, which performs adaptive quadrature.
Numerical root finding methods use iteration, producing a sequence of numbers that hopefully converge towards a limits which is a root. This post only focuses four basic algorithms on root finding, and covers bisection method, fixed point method, Newton-Raphson method, and secant method.
The simplest root finding algorithms is the bisection method. It works when f is a continuous function and it requires previous knowledge of two initial gueeses, u and v, such that f(u) and f(v) have opposite signs. This method is reliable, but converges slowly. For detail, see https://guangchuangyu.github.io/cn/2008/11/bisect-to-solve-equation/ .
Root finding can be reduced to the problem of finding fixed points of the function g(x) = c*f(x) +x, where c is a non-zero constant. It is clearly that f(a) = 0 if and only if g(a) = a. This is the so called fixed point algorithm.
泰勒公式学过微积分都应该知道，可以翻wiki复习一下，https://zh.wikipedia.org/wiki/泰勒公式.
用R简单实现一下：
efv <- function(f, value, variable="x", a=0, eps=0.001) {
#estimate function value using Taylor theorem
assign(eval(variable), a)
fv.old <- eval(f)
k <- 1
repeat {
df <- D(f, variable)
if (df == 0)
break
fv.new <- fv.old + eval(df)*(value-a)^k/factorial(k)
if (fv.new - fv.old < eps)
break
fv.old <- fv.new
f <- df
k <- k+1
}
return (fv.new)
}
翻看了以前写的使用Newton-Raphson Method求一个数的开方，想到其实也可以用中值定理来实现。 中值定理：f(x)是一个连续性的函数，在[u,v]区间内，当c的值位于f(u)和f(v)之间时，至少存在一个点，满足f(x) = c 当f(u)和f(v)一正一负时，那么在[u,v]之间至少有一个根的存在，这个定理本来就是拿来证明根的存在的，但是其实也可以用来求解根。
