上級マクロ経済学

第14回 数値計算の基礎

1 数値計算の基本知識

数値計算

• 経済学や統計学では、関数の最大値（最小値）や連立方程式の解について求めることが非常に多い。
• 最大値や最小値をもたらす値を求めるのが難しい場合が存在する。
  • 関数系が複雑なケース。
• いくつかの連立方程式は、手で解くのが難しかったりする。
  • 方程式が10とか100存在するケース。
  • 関数が陽に解けないケース。

こういった場合には、数値計算を行なって解を求める。

コンピュータができること

• コンピュータは代数の処理ができない。
  • 例えばf(x)=-(x-a)^2, ~(a>0)をグラフとして手描きするなら、なんとなく放物線を書けばいい。
  • f(x)をコンピュータに描いてもらうなら、aに具体的に数値を与えて、次のように順に値を求めていく必要がある。
    • x=-1の時のf(x)の値
    • x=0の時のf(x)の値
    • x=1の時のf(x)の値…

y=-(x-a)^2のグラフ

a=3としたときのグラフを描く。

• コンピュータが関数を描くというのは、「点を繋いでいる」にすぎない。
• 点が多くなれば、本来のなめらかな関数に近づく。
• 最適点の”近似点”を探す、というのが基本。

import numpy as np
import matplotlib.pyplot as plt

## パラメータを設定する。
a = 5

## グラフの点を計算する。
x0 = np.linspace(0,10, 5) # -2から2までを5個に区切ってxをつくる。
x1 = np.linspace(0,10, 10) # -2から2までを10個に区切ってxをつくる。
x2 = np.linspace(0,10, 100) # -2から2までを10個に区切ってxをつくる。
def f(x):
    return -(x-a)**2

plt.plot(x0,f(x0), "-o", label="5 points")
plt.plot(x1,f(x1), "-*", label="10 points")
plt.plot(x2,f(x2), "--", label="100 points")

plt.legend()
plt.show()

アルゴリズム

コンピュータで

• 最大化問題を解く
• 連立方程式の解を求める

ためには、求めたい値を得るための手順に従う。

これはアルゴリズムと呼ばれる。

今回メインで扱うのは、

• 分割法
• ニュートン法

scipyでの最適化

• scipy.optimizeモジュールに、さまざまなアルゴリズムがある。
  • 最小化なら、scipy.optimize.minimize_scalr()。
  • 求根ならsipy.optimize.root_scalar()。
• これらを用いれば、自分で書く必要はない。
• ただし、アルゴリズムを使用するなら、どんなものかはざっくり理解しておくとよい。

2 最大値/最小値

最大化問題

次のような最大化問題を考えよう。

\begin{aligned} \max_{x}&~ f(x) \\ \textrm{s.t.}&~ \underline{x} \le x \le \overline{x} \end{aligned}

• \underline{x} \le x \le \overline{x}の中で、f(x)を最大にするようなxを求める、という問題。

関数の設定

ここからは、 f(x) = -x(x-a)(x-b)(x-c) a = 4,~ b=7,~ c=2 の例で議論する。

最大値をもたらすxの値は解析的に求めることができる。

• f'(x)=0は、xの三次方程式なので、階の公式が存在している。
• ただし、途方もなく複雑な表現になる。
• 数値的に最大値を求めるのが現実的。

総当たり法

一番簡単に思いつくのは、総当たり法（brute-force method）。

• いくつかのxを作る。
  • (x_1, x_2, \cdots, x_n)
  • \underline{x}=x_1,~ \overline{x}=x_n
• f(x_1), f(x_2),\cdotsを計算して、f(x)が最も大きくなる時のx=x_iを探す。
  • np.argmax()関数で、配列の最も大きい値のインデックスを取得する。
• 近似値と真の値とのブレは大きくなる可能性が高い（不安定）。
• 大域的な探索をするならば、大体の値を出せる。

総当たり法のコード

• グリッドを細かく取れば、近似はマシになる。
• 探索領域が小さければ、これでも十分。

import numpy as np
import matplotlib.pyplot as plt

## パラメータを設定する。
a = 4
b = 7
c = 2

## 関数を定義する。
def f(x):
    return -x*(x-a)*(x-b)*(x-c)

## グリッド数 10
x0 = np.linspace(0,8, 10) # 0から8までを10個に区切ってxをつくる。
fx_max0 = np.max(f(x0)) # f(x)の最大の値
ind_max0 = np.argmax(f(x0)) # f(x)の最大値のインデックス
x_max0 = x0[ind_max0] # そのインデックスのxが求めたい値
print("グリッド数10")
print(f"最大のf(x): {f(x_max0)}")
print(f"最大のf(x)をもたらすx: {x_max0}")

## グリッド数 100
x1 = np.linspace(0,8, 100) # 0から8までを100個に区切ってxをつくる。
fx_max1 = np.max(f(x1)) # f(x)の最大の値
ind_max1 = np.argmax(f(x1)) # f(x)の最大値のインデックス
x_max1 = x1[ind_max1] # そのインデックスのxが求めたい値
print("グリッド数100")
print(f"最大のf(x): {f(x_max1)}")
print(f"最大のf(x)をもたらすx: {x_max1}")

plt.plot(x0,f(x0), "-", color="blue", alpha=0.5, label="10 grids")
plt.plot(x_max0, f(x_max0), "*", color="blue", markersize=10)
plt.plot(x1,f(x1), "--", color="green", alpha=0.5, label="100 grids")
plt.plot(x_max1, f(x_max1), "^", color="green", markersize=10)
plt.legend()
plt.show()

グリッド数10
最大のf(x): 45.407712238987976
最大のf(x)をもたらすx: 6.222222222222221
グリッド数100
最大のf(x): 48.088775717483976
最大のf(x)をもたらすx: 5.8989898989899

三分法

• \underline{x} \le x \le \overline{x}単峰性が満たされているならば、三分法を用いることができる。
  • 山が一つだけであればOK。

アルゴリズム:

1. 誤差の許容度\epsilon>0（閾値）を決める。
2. 二箇所にx_0,~x_1~(x_0< x_1)を取る。
   • たいてい、x_0=\underline{x}とx_1=\overline{x}を取ればOK。
3. [x_0,x_1]区間を三等分し、真ん中2つを区切る2点をそれぞれx_{s1}、x_{s2}とする。
4. f(x_{s1})> f(x_{s2})なら、x_1=x_{s2}に置き換える。
5. f(x_{s1})< f(x_{s2})なら、x_0=x_{s1}に置き換える。
6. |x_1-x_0|< \epsilonなら、f(x_0),~f(x_1)のうち大きい法のxの値を解にする。
7. |x_1-x_0|> \epsilonなら、手順3に戻る。

三分法のコード

• 先の関数f(x)で、3\le x \le 8の区間では単峰。
• この区間の最大値を求めよう。

## パラメータを設定する。
a = 4
b = 7
c = 2

## 関数を定義する。
def f(x):
    return -x*(x-a)*(x-b)*(x-c)

def min_tri(func, min, max):
    """
    三分法による最小値探索
    ------------------
    input:
        func : (function) 最小化するxを求めたい関数
        min : (float) xの下限
        max : (float) xの上限

    output:
        x : (float) funcを最小にするx
        fval : (float) fの最小値
    """
    
    x0 = min
    x1 = max
    epsilon = 10**(-6) # 誤差の許容度を決める

    while np.abs(x1-x0) > epsilon:
        
        divide = np.linspace(x0,x1,4) # x0からx1に4つ点をとった配列
        xs1 = divide[1] # 一つ目をxs1
        xs2 = divide[2] # 二つ目をxs2

        if func(xs1) > func(xs2):
            x1 = xs2
        else:
            x0 = xs1
    
    x = np.max([x0,x1])
    return {"x": x, "fval": func(x)}

x = np.linspace(3,8) # xを作る。
res = min_tri(f,3,8) # 最小値を求める
xmax = res["x"] # f(x)を最大にするx
fmax = res["fval"] # 最大のf(x)

plt.plot(x, f(x))
plt.plot(xmax, fmax, "o", markersize=10)

plt.show()

scipyでの最大値探索

• 大抵の統計ソフトは最小値を探索するアルゴリズムを持っている。
  • この場合は、g(x) = -f(x)を作って、g(x)の最小値を求めれば良い。
• pythonで最小値を求めるならば、scipy.optimizeモジュールでminimize_scalar(fun)を用いればよい。
  • ただし、単峰型の関数でなければ、局所的な最小値が求まることに注意が必要。
• ある区間の最小値を求めたいならば、minimize_scalar(fun, bounds=(min, max))を用いる。
  • min\le x \lemaxを求めてくれる。

scipyコード（最大値）

制約なし

import scipy.optimize as opt # optでscipy.optimizeをimport

## パラメータを設定する。
a = 4
b = 7
c = 2

## 最大値がほしい関数を定義する
def f(x):
    return -x*(x-a)*(x-b)*(x-c)
## 最小を返してほしい関数を定義する
def g(x):
    return -f(x)

res = opt.minimize_scalar(g) # 最小値とその時のxを求める
xstar = res.x # gを最小にするx

print(f"x is {xstar}")

x = np.linspace(0,8,100)
plt.plot(x,f(x))
plt.plot(xstar, f(xstar), "o", markersize=10)

x is 0.7763364929761143

制約あり

import scipy.optimize as opt # optでscipy.optimizeをimport

## パラメータを設定する。
a = 4
b = 7
c = 2

## 最大値がほしい関数を定義する
def f(x):
    return -x*(x-a)*(x-b)*(x-c)
## 最小を返してほしい関数を定義する
def g(x):
    return -f(x)

res = opt.minimize_scalar(g, bounds=(3,8)) # 最小値とその時のxを求める（制約付き）
xstar = res.x # gを最小にするx

print(f"x is {xstar}")

x = np.linspace(0,8,100)
plt.plot(x,f(x))
plt.plot(xstar, f(xstar), "o", markersize=10)

x is 5.935362355303199

多変数のケース

• 多変数関数でも、sipy.optimizeのminimize(func, x0)関数を使えば最小値を求めることができる。
• 初期値（x0）が必要なので、適当に与える。 f(x_1,x_2) = x_1 + x_2 - \log (x_1 x_2)

def f(x):
    x1, x2 = x
    return x1 + x2 - np.log(x1*x2)

xmin = opt.minimize(f, x0=[2,2]).x
print(f"fを最小にするx: {xmin}")

x = np.linspace(0.1,5,100)
y = np.linspace(0.1,5,100)
input = np.meshgrid(x,y)

ax = plt.figure().add_subplot(projection = "3d")
ax.plot_surface(input[0],input[1], f(input))
ax.view_init(elev=10, azim=-60)
plt.show()

fを最小にするx: [0.99999996 0.99999996]

3 方程式の求根

求根

• 方程式の解を根という。
• 根を求めることを、求根という。

求根アルゴリズムは、基本は次のような方程式の解を求めるためのアルゴリズム。

f(x) = 0

したがって、方程式を上のような形に整理して、その関数にアルゴリズムを適用する。

例えば、 x^2 = 4 e^{-x} の解となるxを求めるのであれば、 x^2 - 4 e^{-x} = 0 として、 f(x) = x^2 - 4 e^{-x} とする。 このf(x)にアルゴリズムを適用する。 f(x)=0を満たすxが、x^2 = 4 e^{-x}を満たすxになるということがわかる。

f(x)のプロット

def f(x):
    return x**2 - 4*np.exp(x)

x = np.linspace(-2,2,100)
plt.plot(x, f(x))
plt.plot(x, np.zeros(len(x)), color="black", alpha=0.8)
plt.show()

二分法

求根アルゴリズムで最も安定して解を求められるのが、二分法。

1. 誤差の許容度\epsilon > 0（閾値）を決める。
2. 二箇所にx_0, ~ x_1（x_0<x_1）をとる。
   • このとき必ず、f(x_0)の符号とf(x_1)の符号が異なるようにする。
   • f(x)の符号をsign(f(x))とすると、sign(f(x_0)) \neq sign(f(x_1))。
3. [x_0, x_1]区間を二等分し、真ん中の点をx_sとする。
   • x_s = (x_0+x_1)/2でOK。
4. sign(f(x_m))を調べて、
   • sign(f(x_m))=sign(f(x_0))のとき、x_0 = x_mに置き換える。
   • sign(f(x_m))=sign(f(x_1))のとき、x_1 = x_mに置き換える。
5. |x_1 - x_0|>\epsilonならば、手順3に戻る。
6. |x_1 - x_0| \le \epsilonならば、x = (x_0+x_1)/2としてアルゴリズムを終了する。

二分法のコード

def root_bi(func, min, max):
    """
    二分法による求根
    -------------
    input:
        func : (function) 求根したい関数 / f(x)= 0 のf
        min : (float) xの最小値
        max : (float) xの最大値
    output:
        f(x)=0 を満たすx
    """

    epsilon = 10**(-6) # 閾値
    
    if func(min)*func(max) == 0: # 掛け算がゼロならminかmaxのどちらかがf(x)=0の解
        if func(min) == 0:
            return min
        else:
            return max
    elif func(min)*func(max) > 0: # 同符号ならばアルゴリズムが適用できない
        return print("f(min)とf(max)は異符号である必要があります。")
    else: # 異符号なら、アルゴリズムを実行
        x0 = min
        x1 = max
        
        while np.abs(x1-x0) > epsilon:
            xs = (x0+x1)/2 # x0とx1の中間点をxsとする
            
            ## x0かx1をアップデートする。
            if func(xs)*func(x0) > 0: # 
                x0 = xs
            else:
                x1 = xs
        return xs


def f(x):
    return x**2 - 4*np.exp(x)

x_root = root_bi(f, -10,10)
print(f"方程式の解x = {x_root}")
print(f"その時のf(x) = {f(x_root)}")

x = np.linspace(-2, 2, 100)

plt.plot(x,f(x))
plt.plot(x_root, f(x_root), "o", markersize=10)
plt.plot(x, np.zeros(len(x)), color="black", alpha=0.8)
plt.show()

方程式の解x = -1.1342865228652954
その時のf(x) = -2.0603782280304017e-07

勾配法

• 二分法は収束が遅いため、範囲が大きいとそれなりに時間が必要であることが知られている。
  • 後に話す多変数関数（連立方程式）での求根では、求根が困難になる場合もある。
• 勾配法と呼ばれる、関数の傾き（微分の値）を利用した求根がメジャー。
• 最もシンプルなニュートン法を紹介する。

ニュートン法

1. 誤差の許容度\epsilon > 0（閾値）を決める。
2. 初期値x_0を決める。
3. 座標x_0,f(x_0)の接線を描き、接線とx軸との交点を、新たなx_1とする。
4. |f(x_1)-0| = |f(x_1)| \le \epsilonであれば、アルゴリズムを終了する。
5. |f(x_1)| > \epsilonであれば、x_0=x_1と置き換え、手順3に戻る。

接線の求め方

• 関数f(x)の接線の傾きは、その微分係数f'(x)。

• x=x_0において、傾きf'(x_0)で、座標(x_0, f(x_0))を通る直線は、次の式で与えられる。 y = f(x_0) + f'(x_0) (x-x_0)

• この直線がx軸と交わるところは、y=0の時のxの値。したがって、x_1は、 x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}

• f'(x)が手で計算できるなら、それを用いればOK。

• f'(x)が手で計算できない or 面倒な場合は？

  • 数値微分を用いる。

数値微分

微分は、極小区間での関数の傾きを見ているに過ぎない。

• x0での傾きを見たいならば、とても小さい\delta>0を持ってこれば次のような近似ができるはず。 f'(x_0) \equiv \frac{f(x_0+\delta) + f(x_0)}{\delta}

• ただし、近似を行うという性質上、数値微分では次のように中心化した微分が用いられる。 f'(x_0) \equiv \frac{f(x_0+\delta) + f(x_0-\delta)}{2\delta}

• 曲がりのきつい関数であれば近似制度は低くなる。

  • \deltaをどれほど小さく取るか、という問題にも関わる。

数値微分比較

• pythonでは、scipy.derivativeモジュールのderivative()が用意されている。

import scipy.differentiate as spdif

def f(x):
    return np.exp(x)

def deriv(func, x0, delta=10**(-6)):
    """
    数値微分アルゴリズム
    -----------------
    input:
        func : (function) 微分したい関数
        x0 : (float) funcを微分する位置
        delta : (float/ default=10**(-6)) 微分に用いる微小区間
    output:
        (float) funcをx0で微分した値
    """
    return (func(x0+delta) - func(x0-delta))/(2*delta)


print(f"e^xを微分。真の値は{np.exp(1)}")
print(f"自分で定義 : {deriv(f,1)}")
print(f"scipy.differentiate.derivative : {spdif.derivative(f, 1).df}")

e^xを微分。真の値は2.718281828459045
自分で定義 : 2.718281828295588
scipy.differentiate.derivative : 2.7182818284590127

ニュートン法のコード

• 数値微分を用いてニュートン法を実装する。

## 数値微分のアルゴリズムを定義する
def deriv(func, x0, delta=10**(-6)):
    """
    数値微分アルゴリズム
    -----------------
    input:
        func : (function) 微分したい関数
        x0 : (float) funcを微分する位置
        delta : (float/ default=10**(-6)) 微分に用いる微小区間
    output:
        (float) funcをx0で微分した値
    """
    return (func(x0+delta) - func(x0-delta))/(2*delta)

## ニュートン法のアルゴリズム
def root_newton(func, x0):
    """
    ニュートン法による求根
    -------------
    input:
        func : (function) 求根したい関数 / f(x)= 0 のf
        x0 : (float) xの初期値
    output:
        f(x)=0 を満たすx
    """

    epsilon = 10**(-6) # 閾値
    
    while np.abs(func(x0)) > epsilon:
        
        fp = deriv(func, x0)
        x1 = x0 - func(x0)/fp
        x0 = x1
        
    return x0


def f(x):
    return x**2 - 4*np.exp(x)

x_root = root_newton(f, 0)
print(f"方程式の解x = {x_root}")
print(f"その時のf(x) = {f(x_root)}")

x = np.linspace(-2, 2, 100)

plt.plot(x,f(x))
plt.plot(x_root, f(x_root), "o", markersize=10)
plt.plot(x, np.zeros(len(x)), color="black", alpha=0.8)
plt.show()

方程式の解x = -1.1342868185419477
その時のf(x) = 8.45145682992765e-07

scipy.optimizeでの求根

• scipy.optimizeモジュールのroot_scalr()を用いれば、求根は簡単にできる。
  • root_scalr(func,bracket=(a,b))とすると、二分法で求根してくれる。
  • root_scalr(func,x0=a)とすると、初期値aの勾配法で求根してくれる。
    • しかも、ニュートン法よりもアップデートされたアルゴリズムを使ってくれる。
    • 数値微分なども勝手にやってくれる。
  • 返り値は色々な値をセットにしたものなので、そのうちのroot成分を取る（以下のコード参照）。

import scipy.optimize as opt

def f(x):
    return x**2 - 4*np.exp(x)

## 二分法
res_bi = opt.root_scalar(f, bracket=(-10,10))
x_root_bi = res_bi.root

## 勾配法
res_grad = opt.root_scalar(f, x0=0)
x_root_grad = res_grad.root

print(f"二分法の解x = {x_root_bi}")
print(f"勾配法の解x = {x_root_grad}")

二分法の解x = -1.1342865808195677
勾配法の解x = -1.1342865808195677

連立方程式の解

• 関数の変数がn変数あるときに、n本の連立方程式があれば、解が定められる。

例： \begin{aligned} f(x_0, x_1) &= x_0 + 0.5(x_0 - x_1)^3 - 1 = 0, \\ g(x_0, x_1) &= 0.5(x_1 - x_0)^3 + x_1 = 0 \end{aligned}

3次元

2次元(z=0)

[0.8411639 0.1588361]

連立方程式の求根

• このような連立方程式についても、sicpy.optimizeのroot(func, x0)関数を使えば解が得られる。
  • たいてい、勾配法を使う。
• funcの形は、次のような形になる。

input:
    x : (array) 2次元の配列。
output:
    (list) [f(x[0],x[1]), g(x[0],x[1])]

• 例を見てみよう。

def fun(x):
    return [x[0]  + 0.5 * (x[0] - x[1])**3 - 1.0,
            0.5 * (x[1] - x[0])**3 + x[1]]

pit = opt.root(fun, x0=[0,0]).x
print(f"[x0, x1] = {pit}")

[x0, x1] = [0.8411639 0.1588361]

まとめ

• 最小値を求める、方程式の解を求める、の二つが数値計算の基礎。
• それなりにテクニカルだが、基本的なアルゴリズムを覚えたら、あとはそれらの改良版が出ているんだな、という認識で良い。
• 我々はエンドユーザーなので、用意されているパッケージがあればそれを用いるのが良い。
  • 数値微分についてなどは別。
• 次週はソローモデルを紹介して、数値解を求める。