文系の学生などで特に目的無くプログラミングを覚えたい人には、Rが最適です。統計解析に使わなくても、ちょっと計算した気になるプログラミングの練習をしてみましょう。
1. よくある近似解を求める
関数f(x)=0があるとして、テーラー展開をしてみます。f(x)=f(x0)+f'(x0)(x-x0)+o(x-x0)=0になり、o(x-x0)は十分小さいと無視すると、x=x0-f(x0)/f'(x0)と整理できます（Wikipedia）。ここでx2+x-5=0のときのxをRを使って求めてみましょう。
# 方程式を設定する
f <- expression(x^2+x-5)
# 一階微分を求める
df <- D(f, "x")
# 初期値を定義
x <- 10
# 近似計算を実行
x <- x - eval(f)/eval(df)
print(paste("x:",x))とすると、x=5.00が得られます。まだ答えには遠い。しかし、x <- x - eval(f)/eval(df) を繰り返すと2.72...となり、続けて5回も繰り返すと1.79...で動かなくなります。この時点でx^2+x-5を計算すると0になるので、xの解になっている事が分かります。この漸近的な近似解を求めるアルゴリズムを、ニュートン・ラフソン法(Newton-Raphson method)と言いいます。
2. ループで手抜きをする
ニュートン・ラフソン法がプログラミング言語の紹介として良い特性があるのは、何回も繰り返す操作があるからです。ちょっと面倒なので、次のようなループで処理を行いましょう。
f <- expression(x^2+x-5);
df <- D(f, "x");
x=10;
# ループは10回
for(i in 1:10){
x <- x - eval(f)/eval(df);
print(sprintf("%02d times: x=%2.5f",i,x));
}
2.1.補足
目的関数の最小値を計算するnlm()関数を使えば、もっと簡潔に実行できます。
nlm(function(p){
(p^2+p-5)^2
}, 10)
なお、初期値10と-10で演算結果が異なりますが、どちらも正しいです。
3. ニュートン・ラフソン法で連立方程式を解く
経済学に限らず多変量のモデルを扱う事が多く、一次元のモデルでは実用性が低いように感じられるかも知れません。そこで二次元に応用してみましょう。数学的には行列の微分になるだけで、テイラー展開がベースになっているのは変わりませんが、Rのコードはヤコビアンを作ったり逆行列を求めたりするので、少し複雑になります。
# 連立方程式を設定
f1 <- expression(x^2 + x - y - 4)
f2 <- expression(x - 2*y)# xとyで一階微分を作っておく
g1 <- deriv(f1, c("x", "y"))
g2 <- deriv(f2, c("x", "y"))# 初期値を設定
x <- 110
y <- -10for(i in 1:10){
# (x,y)を2x1行列にする
m <- matrix(c(x, y), 2, 1)
# (x,y)で評価したf1とf2の値を2x1行列にする
f <- matrix(c(eval(f1), eval(f2)), 2, 1)
# ヤコビアンを作る
# グラディエント・ベクトルg1とg2から行列形式で値を得るにはattr()が必要なので注意
j <- rbind(attr(eval(g1), "gradient"), attr(eval(g2), "gradient"))
# 行列mから、ヤコビアンの逆行列と評価式を乗じたものを引く
m <- m - solve(j)%*%f
# 行列mを(x, y)に展開しておく
x <- m[1]; y <- m[2]
print(sprintf("%d times: (x,y)=(%f,%f)", i, x, y))
}
9回目ぐらいの計算で、x=1.765564、y=0.882782が得られると思います。三次、四次と拡張していけば、ちょっと手計算が難しいモデルでも、数値解を得ることができるかも知れません。
何をやっているのかイメージが沸かない人のために、初期値を(15,15)にした場合のニュートン・ラフソン法による近似を4順目までプロットしてみました。初期値から交点へ接近していく過程が分かると思います。
なおループを繰り返しても交点に接近しなくなる収束点に辿り着くのは7順目です。
3.1. 補足
nleqslv()関数で簡単に計算ができます。以下は利用例。
library("nleqslv")
f <- function(p){ return(c(p[1]^2 + p[1] - p[2] - 4, p[1] - 2*p[2])) }
nleqslv(c(15, 15), f)
