ウォリスの公式

PythonでWallisの公式を書いてみました。Wallisの公式は円周率を多項式を利用して近似するものです。元ネタはPython Scientific lecture noteの英語版(2012年8月リリース)19ページの練習問題になります。実行結果を見てみると、引数の値を大きくするにつれて、円周率に近づいている様子が確認できます。ちなみにnumは数字というよりもnumerator(分子)のつもりです、念のため。

def wallis_formula(n):
    mypi = 1
    for i in range(1, n):
        num = 4*(i**2)
        mypi *= float(num)/(num-1)
    return 2 * mypi

'''
>>> wallis_formula(1000)
3.1408069608284657
>>> wallis_formula(10000)
3.1415141108281714
>>> wallis_formula(100000)
3.141584799578707
>>> wallis_formula(1000000)
3.1415918681913633
'''

ニュートン法による極値計算

テキスト「これなら分かる最適化数学」の86ページにあるニュートン法をPythonで実装した。通常のニュートン法は関数f(x)に対し、f(x)=0なるxを求める反復計算を行うものであるが、反復の式を少し変更するだけで、極値(f'(x)=0なるx)を求めることができる。詳しい理論はテキストに譲るが、x <- x - f'(x)/f''(x)とし、収束するまで繰り返せば、極値にたどり着くことができる。

ニュートン法で極値を求めるにあたり、注意しなければならないのは、コードの更新式からも見当がつくが、途中で2回微分が0にならないようにすることである。

以下、f(x)=x^3-2*x^2+x+3に対して、極値が求まるか試してみた。なお、f(x)が極値をとるxの値は解析的に解くこともでき、x=1、1/3の2つとなる。

# -*- encoding: utf-8 -*-

'''
ニュートン法による極値の計算
2回微分が存在する場合、勾配法よりも効率よく極値を計算できる
'''

from sympy import * # 微分操作のために必要
import itertools as it # mapで複数の引数を与えたいので利用する

def newton(init, vari, func, delta=10**-6):
    '''
    funcは極値を探す対象の関数、variは微分する変数、deltaは収束判定のための定数
    返す値は見つかった極値の1つである
    '''
    df = diff(func, vari)
    d2f = diff(df, vari)
    x = init
    while True:
        xvar = x
        # 2回微分の微分係数が0にならないように注意が必要
        x = xvar - float(limit(df, vari, xvar)) / float(limit(d2f, vari, xvar))
        if abs(x - xvar) < delta:
            break
    return round(x, 5) # deltaの値によって有効桁数が変わるが、簡単のため5で決めうち

x = Symbol('x')
f = x**3 - 2*x**2 + x + 3
start_num = range(-5, 5) # 初期値を複数設定しておく
print map(newton, start_num, it.repeat(x, len(start_num)),
          it.repeat(f, len(start_num)))

最後のmapを利用することで、一度に複数の結果を実行、表示することができるようになる。mapに与える引数については、繰り返し実行できる形式にしておく必要がある。

実行結果は、[0.33333, 0.33333, 0.33333, 0.33333, 0.33333, 0.33333, 1.0, 1.0, 1.0, 1.0]となる。それぞれ初期値を-5、-4、...、4に変更しながら実行した結果である。すべて1/3と1のうち、近いほうの極値に収束していることが確認できる。

勾配法を用いた最適化計算

テキスト「これなら分かる最適化数学」の81ページにある勾配法(gradient method)アルゴリズムをpythonで実装した。勾配法とは、関数f(x)を最大にするxの値に近いと思われる値を初期値とし、必ず関数値f(x)が増加するようにし、かつステップ幅をなるべく大きくなるようその幅を調整することにより、できるだけ高速に関数を最大にするxを求めるというものである。

なお、今回扱う勾配法においては、変数の値がとりうる範囲について、関数の微分係数が0になるのは1箇所であると仮定する(複数極値があると、どこに向かうか分からないため)。また、微分できない箇所がある関数では不都合が生じるので、これも考えないことにする(例えばf(x)=|x|のx=0の点のように、尖っている状況下などがある)。

関数の最大値を求めるにあたっては、微分係数が0になるという性質を利用している。このため、微分という操作を行う必要であるが、この問題はsympyモジュールを利用することによって解決できる。sympyモジュールはこちらからダウンロードできる。なお、記事執筆時点のsympyのバージョンは0.7.1である。

以下、「2次関数f(x)」と「正規分布の確率変数と関係がある部分のみ取り出した関数g(x)」の2つの例を用いて、勾配法で正しく関数を最大化するようなxを推定できるか試してみた。結果としては、(制約が強いこともあるが)ほぼ正確に推定できているということになった。

# -*- encoding: utf-8 -*-

from sympy import *

# 符号関数
def sign(value):
    if value > 0: return 1
    if value < 0: return -1
    return 0

# 勾配法を計算するメインルーチン
# initは初期推定値、funcは最大値を求めたい関数、stepは初期移動量、
# epsは収束した(ほぼf'(x)が0となった)と判定する閾値である。
# 返り値は最大値となるxの値に相当する。
def gradient_method(init, func, step = 10**-7, eps = 10**-6):
    # diffは微分の関数、これで第一引数の関数をxで微分したことになる
    # 事前にxをSymbol関数でシンボル化しておく必要がある
    df = diff(func, x)
    xx, h = init, step
    while True:
        h = sign(limit(df, x, xx)) * abs(h)
        X, Xdash = xx, xx + h
        # 関数に値を入れるという方法が見つからないので
        # 極限の値を求めている
        if limit(func, x, X) < limit(func, x, Xdash):
            while limit(func, x, X) < limit(func, x, Xdash):
                h *= 2.0
                X, Xdash = Xdash, X + h
            xx, h = X, h / 2.0
        else:
            while limit(func, x, X) > limit(func, x, Xdash):
                h /= 2.0
                Xdash -= h
            xx, h = Xdash, 2.0 * h
        if abs(limit(df, x, xx)) <= eps:
            break
    return xx

x = Symbol('x') # xを記号として扱う
f = -x**2 + 2*x + 3
g = exp(-(x-2)**2)
print gradient_method(0, f) # 1.0000002、なお正解は1である
print gradient_method(0.5, g) # 2.0000002、なお正解は2である

バームクーヘン型積分

よくある1変数積分の問題に、ある関数で囲まれる面をy軸の周りに1回転させてできる体積を求めよというものがあります。そのような問題を解く際に用いられる手法としてバームクーヘン型積分というものが知られています。バームクーヘン(ドイツ語で木のケーキという意味です)の形状は軸の周りに線が入っているのが特徴です。バームクーヘン型積分というのは、その線を一つ一つかき集めて体積とみなす計算方法です。

具体的にy=f(x) (0<=x<=a)という形が与えられていると、この曲線とx軸で囲まれた面をy軸回りで1回転させたときの体積は次式のようになります。

2*PI*x*yというのは、底面の半径がx、高さがf(x)の円柱の側面積に相当します。側面積をxが0からaまでの範囲で集まると先述のように体積となるわけです。

余談ですが、この問題は東大が大学受験の問題として出したのがきっかけで広まったそうです。