8. SciPy¶

8.2.1. Random Variables and Distributions¶

この小節では，ベータ分布を使ってみます．まず， numy.random.beta でベータ分布に従うランダムな数字をだしてみましょう．

>>> import numpy as np
>>> np.random.beta(5, 5, size=3) #ベータ分布
array([ 0.5406392 ,  0.25673906,  0.41268416])

ところで，ベータ分布とは，以下の様な確率密度関数に従うような分布でした．

この関数の分母にあるのは，ベータ関数ですね．

独立に一様分布 \(U(0,1)\) に従う \(p+q-1\) 個の確率変数を大きさの順に並べ替えたとき，小さい方から p 番め（大きい方からは q 番目）の確率変数 X の分布がベータ分布 \(B(p,q)\) となります。

以上のような，形でベータ分布に従う数字を生成することが出来ますが，特定のパラメータと特定の分布に従うような数字を生成する必要がある場合もあるとおもいます．

例えば，以下のような，図を出したいときです．

黒い実線の曲線はベータ分布の確率密度関数です(理論値)．青い棒グラフのようなものは，先ほどのようにして，ランダムにベータ分布から数字を生成して，プロットしていったものです(観測値)．黒い実線が補助線になっていて，青い観測値がベータ分布にしてがっているような雰囲気を感じることができます．

このような図を生成するのに用いたコードは以下のようなものです．

import numpy as np
from scipy.stats import beta
from matplotlib.pyplot import hist, plot, show
%matplotlib inline
q = beta(5, 5)      # Beta(a, b), with a = b = 5 scipyじゃないと出ない？？？
obs = q.rvs(2000)   # 2000 observations
hist(obs, bins=40, normed=True)
grid = np.linspace(0.01, 0.99, 100)
plot(grid, q.pdf(grid), 'k-', linewidth=2)
show()

このコードの中では， q =beta(5, 5) で観測値を作る仮定で，いわゆる， rv_frozen object というものが作られています．

この frozen を使うことで，分布を固定して，観測値を生成することが出来ます．

この他にも， q = beta(5,5) にたいして， q.cdf() , q.pdf() , q.ppf() , q.mean() などが使うことが出来ます．

8.2.2. Other Goodies in scipy.stats¶

今まで説明した以外にも， SciPy には有益な 機能があります．

例えば， scipy.stat.linregress では，線形回帰分析を行うことができます．

>>> import numpy as np
>>> from scipy.stats import linregress

>>> x = np.random.randn(200)
>>> y = 2*x + 0.1*np.random.randn(200)
>>> gradient, intercept, r_value, p_value, std_err = linregress(x, y)
>>> gradient, intercept
(2.0034904702816632, -0.0038496543788882898)

と言った具合です．

8.3.1. Bisection¶

Bisection（二分法）とは，以下の様な簡単なアルゴリズムで根を見つける方法です．このアルゴリズムを理解するために，以下の様な簡単なゲームを考えてみましょう．

プレイヤーAが，プレイヤーBが1から100の間から選んだ数字を当てるゲームをします．この時，以下の様なアルゴリズムでAが数字を特定していったとき，このアルゴリズムをBisectionと呼びます．

• 選んだ数字は50以下か？

  • もし50以下なら，それは25以下か
  • もし75以下なら，それは75以下か

これが，Bisectionです．つまり，解が存在する定義域の範囲を半分にしていって，解を特定していく方法です．この，アルゴリズムは連続な増加関数に対して有効に働きます．

以下に，このアルゴリズムを実装した関数のコードを掲載します．

def bisect(f, a, b, tol=10e-5):
"""
Implements the bisection root finding algorithm, assuming that f is a
real-valued function on [a, b] satisfying f(a) < 0 < f(b).
"""
lower, upper = a, b

while upper - lower > tol:
        middle = 0.5 * (upper + lower)
        # === if root is between lower and middle === #
        if f(middle) > 0:
        lower, upper = lower, middle
        # === if root is between middle and upper  === #
        else:
        lower, upper = middle, uppera

return 0.5 * (upper + lower)

ちなみに，関数の引数の中に入っている， tol=10e-5 は許容する誤差の範囲です．Bisectionを繰り返していって，どの程度まで絞り込めたら，根の探索をやめるかをここで指定しています．

もちろん，PythonにはこのBisectionの関数が備わっています．

>>> from scipy.optimize import bisect
>>> f = lambda x: np.sin(4 * (x - 0.25)) + x + x**20 - 1
>>> bisect(f, 0, 1)
0.4082935042797544

8.3.2. The Newton-Raphson Method¶

次に説明するのは，The Newton-Raphson Method（一般的に，ニュートン法と呼ばれているようです）です．この，手法も頻繁に用いられる方法です． SciPy では scipy.newton として用います．Bisectionとは異なり，The Newton-Raphson Methodは，根を探す関数の傾きを使って根を求めます．

The Newton-Raphson Methodは諸刃の剣のような性質を持っています．The Newton-Raphson Methodは，関数が根を探す範囲で微分可能であれば，Bisectionよりも早く根を見つける事ができます．しかし，例えば，不連続な点や，微分不可能で傾きが求められない場合，根を見つける事ができません．

では，このThe Newton-Raphson Methodをつかって，Bisectionでやった時と同じ関数の根を見つけてみましょう．

>>> from scipy.optimize import newton
>>> newton(f, 0.2)
0.408293504279
>>> newton(f, 0.7)
0.70017

8.3.3. Hybrid Methods¶

前小節でみたように，The Newton-Ranpson methodは速く根を求める事ができますが，用いることが出来ない場合もあります(non-robast)．一方で，Bisectionは，比較的適用できる関数が多くあります(robast)．

なので，実際の場面では，このどちらかの２つを状況に合わせて使い分けるのがいいでしょう．つまり，最初に速い手法を試して，その結果を確認し，それが正しくないと判断できたら，よりロバストな手法を適用すればいいわけです.

実際, scipy.optimize の brentq ではこのようなhybrid-methodが用いられています．

>>> from scipy.optimize import brentq
>>> print brentq(f, 0, 1)
0.408293504279

8.3.4. Fixed Point¶

Scipyを用いて，不動点を探すことも出来ます．

試してみましょう，

>>> from scipy.optimize import fixed_point
>>> fixed_point(lambda x: x**2, 10.0)  # 10.0 is an initial guess
1.0

8.3.5. Optimization¶

Scipyを用いて最小値を求める事もできます．同様に，最大値も求めることが出来ます．なぜなら，ある関数fの最小値が求まるということは，-fの最小値も求めることもでき，その値はfの最大値でもあるからです．

この最小化は，根の探索とよく似ています．つまり，滑らかな関数であれば，その関数の最適な値を求めるということはその関数の１階微分の根を探すことと同じだからです．

特に，目的関数に関する情報を持ってないのであれば，hybrid-methodを用いるのが良いでしょう.スカラー関数の，制約付き最適化問題の場合，hybrid-methodである fminbound が便利です．

>>> from scipy.optimize import fminbound
>>> fminbound(lambda x: x**2, -1, 2)  # Search in [-1, 2]
0.0