SciPy http://www.scipy.org/ は様々な科学技術計算のための実装を含むパッケージ（モジュール群）である．
- optimize: 非線形最適化と方程式の根
- spatial: 計算幾何学
- stats: 確率・統計
- interpolate: 補間
- integrate: 積分
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return x / (1 + x**2)
x = np.linspace(-10, 10, 100)
y = f(x)
plt.plot(x, y)
from scipy.optimize import minimize_scalar
print("Brent=========================")
print(minimize_scalar(f, method="Brent"))
print("Golden=========================")
print(minimize_scalar(f, method="Golden"))
print("Bounded=========================")
print(minimize_scalar(f, bounds=(-2.0, -1.5), method="Bounded"))
SciPyのoptimizeサブモジュールの minimize(f,x0)
複数の変数をもつ関数fを初期解x0からの探索で最小化
引数methodで探索のためのアルゴリズムを設定
'Nelder-Mead'：Nelder–Mead法(単体法)
'Powell'：Powellの共役方向法
'CG'：共役勾配法(conjugate gradient)
'BFGS'：Broyden–Fletcher–Goldfarb–Shanno(BFGS)法
'Newton-CG'：Newton共役勾配法
'L-BFGS-B'：記憶制限付きBFGS法
'TNC'：打ち切りNewton共約勾配法
'COBYLA'：Constrained Optimization BY Linear Approximation
'SLSQP'：Sequential Least SQuares Programming
'dogleg'：ドッグレッグ信頼領域法
'trust-ncg'：信頼領域Newton共約勾配法
他の引数
jac : 勾配ベクトル(Jacobian)
hess : Hesse行列(Hessian)
bounds :上下限（限界値）制約
constraints: 一般の制約式（$=0$ もしくは $\geq 0$）
制約の種類は文字列 'type' をキーとし，値を制約の種類を表す文字列 ('eq', 'ineq')．
制約関数は文字列 'fun'をキーとし，制約を表す関数を値
Jacobianは文字列 'jac'をキーとし，Jacobianを表す関数を値
制約の例：不等式 $x_0 +2x_1 \geq 10$ を表す制約
{'type': 'ineq', 'fun': lambda x: x[0] + 2*x[1] -10}
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import LogNorm
import matplotlib.cm as cm
ax = Axes3D(plt.figure())
s = 0.05
X = np.arange(-2, 2.0 + s, s)
Y = np.arange(-1, 3.0 + s, s)
X, Y = np.meshgrid(X, Y)
Z = (1.0 - X) ** 2 + 100.0 * (Y - X * X) ** 2
ax.plot_surface(
X, Y, Z, rstride=1, cstride=1, norm=LogNorm(), cmap=cm.jet, linewidth=0.2
)
None
import plotly.graph_objs as go
import plotly
X, Y = np.mgrid[-2:2:0.05, -1:3:0.05]
Z = (1.0 - X) ** 2 + 100.0 * (Y - X * X) ** 2
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=np.log(Z + 0.01)))
# fig.show()
# plotly.offline.plot(fig); #Jupyter Labの場合
from scipy.optimize import minimize, rosen, rosen_der, rosen_hess
x0 = [2.0 for i in range(5)]
res = minimize(rosen, x0, method="Nelder-Mead")
print(res.x)
res = minimize(rosen, x0, jac=rosen_der, method="BFGS")
print(res.x)
res = minimize(rosen, x0, jac=rosen_der, hess=rosen_hess, method="trust-ncg")
print(res.x)
ax = Axes3D(plt.figure())
c = {1: 1.5, 2: 2.25, 3: 2.625}
s = 0.5
X = np.arange(-5, 5.0 + s, s)
Y = np.arange(-5, 5.0 + s, s)
X, Y = np.meshgrid(X, Y)
Z = sum((c[i] - X * (1 - Y**i)) ** 2 for i in range(1, 4))
ax.plot_surface(
X, Y, Z, rstride=1, cstride=1, norm=LogNorm(), cmap=cm.jet, linewidth=0.2
)
c = {1: 1.5, 2: 2.25, 3: 2.625}
X, Y = np.mgrid[-5:5:0.05, -5:5:0.05]
Z = sum((c[i] - X * (1 - Y**i)) ** 2 for i in range(1, 4))
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=np.log(Z + 0.01)))
# fig.show();
# plotly.offline.plot(fig); #Jupyter Labの場合
ax = Axes3D(plt.figure())
s = 0.1
X = np.arange(-3, 3.0 + s, s)
Y = np.arange(-3, 3.0 + s, s)
X, Y = np.meshgrid(X, Y)
Z = (
-20 * np.exp(-0.2 * np.sqrt(0.5 * (X**2 + Y**2)))
- np.exp(0.5 * (np.cos(2 * np.pi * X) + np.cos(2 * np.pi * Y)))
+ np.exp(1)
+ 20
)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet, linewidth=0.2)
None
X, Y = np.mgrid[-3:3:0.05, -3:3:0.05]
Z = (
-20 * np.exp(-0.2 * np.sqrt(0.5 * (X**2 + Y**2)))
- np.exp(0.5 * (np.cos(2 * np.pi * X) + np.cos(2 * np.pi * Y)))
+ np.exp(1)
+ 20
)
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=Z))
# fig.show()
# plotly.offline.plot(fig); #Jupyter Labの場合
X, Y = np.mgrid[-3:3:0.05, -3:3:0.05]
Z = X**2 + X * Y + Y**2
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=Z))
# fig.show()
# plotly.offline.plot(fig); #Jupyter Labの場合
X, Y = np.mgrid[-3:3:0.05, -3:3:0.05]
Z = np.sqrt(1 + X**2 + Y**2)
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=Z))
# fig.show()
# plotly.offline.plot(fig); #Jupyter Labの場合
X, Y = np.mgrid[-3:3:0.05, -3:3:0.05]
Z = X**3 - 3 * X * Y + Y**3
fig = go.Figure(go.Surface(x=X, y=Y, z=Z, surfacecolor=Z))
# fig.show()
# plotly.offline.plot(fig); #Jupyter Labの場合
制約付き最適化の例
半径 $\sqrt{2}$ の円内の制約の下で $x+y$ を最大化
COBYLA (線形近似法;Constrained Optimization BY Linear Approximation）で求解
SLSQP（逐次最小2乗法;Sequential Least SQuares Programming）で求解
目的関数や制約式は，普通の関数(defで定義）で書いても，lambda関数として書いても良い．以下の例では2通りの書き方で目的関数を書いている．（SciPyでの目的関数は最小化なので，-1を乗じてある．）
制約は辞書を用いて定義する． 辞書には，以下のキーと値の対を入れる．
- 'type': 制約のタイプで等式のときは 'eq'，不等式のときには 'ineq'を入れる．
- 'fun': 制約を表す関数を入れる． 制約は必ず，$0$ 以上の制約に直してから入れる．
以下の例では， $x^2 + y^2 \leq 2$ を $-x^2 -y^2 +2 \geq 0$ に直してから入れている．
その後に，minimmize関数で最小化するが，制約が複数ある場合には，制約（辞書）のリストを constraints引数に入れる．
from scipy.optimize import minimize
def f(x):
return -x[0] - x[1]
bnds = ((0, None), (0, None)) # 上下限制約の入れ方
f = lambda x: -x[0] - x[1]
con = {"type": "ineq", "fun": lambda x: -x[0] ** 2 - x[1] ** 2 + 2}
res = minimize(fun=f, x0=[2, 2], method="COBYLA", constraints=con)
print("COBYLA========= \n", res)
res = minimize(fun=f, x0=[2, 2], method="SLSQP", constraints=con)
print("SLSQP========== \n", res)
問題（制約付きWeber問題）
以下のような位置と人数をもった7件の家がある．
| 家 | x座標 | y座標 | 人数 |
|---|---|---|---|
| 1 | 24 | 54 | 2 |
| 2 | 60 | 63 | 1 |
| 3 | 1 | 84 | 2 |
| 4 | 23 | 100 | 3 |
| 5 | 84 | 48 | 4 |
| 6 | 15 | 64 | 5 |
| 7 | 52 | 74 | 4 |
みんなで出資して新しく井戸を掘ることになったが， 話し合いの末 「各家が水を運ぶ距離 $\times$ 各家の水の消費量」の総和が最小になる場所に 井戸を掘ることにした． ただし, 各家の水の消費量は人数に比例するものとする．
どこを掘っても水が出るものとしたとき，どのようにして掘る場所を決めれば良いだろうか？
$(50,50)$ から半径 $40$ の円状の領域が砂漠になっていて，そこには井戸が掘れないときにはどうしたら良いだろうか？
家の集合を $H$，家 $i (\in H)$ の位置を$(x_i, y_i)$ とする． 井戸の位置を $(X, Y)$ とすれば，家 $i$ から井戸までの距離は $$ \sqrt{(x_i - X)^2 + (y_i - Y)^2} $$ である．家 $i$ が $1$日に必要とする水の量（人数）を $w_i$ としたとき， $$ \sum_{i \in H} w_i \sqrt{(x_i-X)^2 + (y_i-Y)^2} $$ を最小にする $(X, Y)$ を求めよ．
定式化： $$ \begin{array}{ll} minimize & \sum_{i \in H} w_i \sqrt{(x_i-X)^2 + (y_i-Y)^2} \\ s.t. & (50-X)^2 +(50-Y)^2 \geq 40^2 \end{array} $$
x = [24, 60, 1, 23, 84, 15, 62]
y = [54, 63, 84, 100, 48, 64, 74]
w = [2, 1, 2, 3, 4, 5, 4]
問題（ポートフォリオ最適化）
100万円のお金を5つの株に分散投資したいと考えている． 株の価格は，現在はすべて $1$株あたり $1$円だが， 証券アナリストの報告によると，それらの株の $1$年後の価格と標準偏差はそれぞれ 以下の表のように確率的に変動すると予測されている．
| 株式 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| 期待値 ($r_i$) | 1.01 | 1.05 | 1.08 | 1.10 | 1.20 |
| 標準偏差 ($\sigma$) | 0.07 | 0.09 | 0.1 | 0.2 | 0.3 |
目的は$1$年後の資産価値を最大化することである． しかしながら，よく知られているように，1つの株式に集中投資するのは 危険であり，大損をすることがある． 期待利回りを $\alpha$%以上としたとき，標準偏差の2乗和を最小にするように投資するにはどうすれば良いだろうか？
ポートフォリオ最適化（定式化）:
$$ \begin{array}{lrcl} minimize & \sum_{i=1}^n \sigma_i^2 x_i^2 \\ s.t. & \sum_{i=1}^n r_i x_i &\geq& \alpha \\ & \sum_{i=1}^n x_i &=& 1 \\ & x_i &\geq& 0 \ \ \ \forall i=1, 2, \ldots, n \end{array} $$問題（学区割り当て問題） （難）
- 2 つの学校（定員100人）を作る
- 学生たちの歩く距離の合計を最小
- 新しい学校の位置と割当
| 学区 | x座標 $x_i$ | y座標 $y_i$ | 人数 $n_i$ |
|---|---|---|---|
| 1 | 0 | 0 | 40 |
| 2 | 0 | 100 | 40 |
| 3 | 100 | 0 | 40 |
| 4 | 100 | 100 | 40 |
| 5 | 50 | 50 | 40 |
学区割り当て問題(定式化）
$$ \begin{array}{l l l} minimize & \sum_{i \in I} \sum_{j=0}^1 z_{ij} \sqrt{ (x_i-X_j)^2+(y_i-Y_j)^2 } & \\ s.t. & \sum_{i \in I} z_{ij} \leq 100 & \forall j=0,1 \\ & \sum_{j=0}^1 z_{ij} = n_i & \forall i \in I \\ & z_{ij} \geq 0 & \forall i \in I; j=0,1 \end{array} $$- 非凸関数の最適化問題
- 適切な初期解を与えないと妥当な解が得られない可能性
問題（複数品目の経済発注量問題）
古典的な経済発注量モデルの複数品目への拡張（容量制約付き）
- $I$ : 品目の集合
- $w_i$ : 品目 $i~(\in I)$ の大きさ
- $W$ : 倉庫に置ける量の上限．
- $F_i$ : 品目 $i$ の発注費用
- $h_i$ : 品目 $i$ の在庫保管費用
- $d_i$ : 品目 $i$ 需要量
容量制約を破らない条件の下で，総費用が最小になる発注方策を求める．
複数品目の経済発注量問題（定式化）
$T_i$：品目 $i$ の発注間隔（発注頻度の逆数）
$$ \begin{array}{l l l} minimize &\sum_{i \in I} \frac{F_i}{T_i} + \frac{h_i d_i T_i}{2} & \\ s.t. & \sum_{i \in I} w_i d_i T_i \leq W & \\ & T_i > 0 & \forall i \in I \end{array} $$数値は自分で適当に定めて解け．
from scipy.optimize import curve_fit
def f(x, a, b, c):
return a * x**2 + b * x + c
xdata = np.linspace(-10, 10, 20)
ydata = f(xdata, 1, 1, 1) + np.random.normal(size=len(xdata))
param, cov = curve_fit(f, xdata, ydata)
print(param)
plt.plot(xdata, ydata, "ro")
plt.plot(xdata, f(xdata, param[0], param[1], param[2]))
None
def f(x):
return x**2 - 4 * np.sin(x)
x = np.linspace(-1, 3, 100)
plt.plot(x, np.linspace(0, 0, 100)) # x=0
y = f(x)
plt.plot(x, y)
None
from scipy.optimize import root, brentq
print("Brentq=", brentq(f, 1.0, 3.0))
print("Root finding ======\n ", root(f, 2.5))
$K$-d木
KDTreeは$K$次元の点データに対する$K$-d木とよばれるデータ構造のクラスである． ここで$K$-d木($K$-dimensional tree)とは， $K$次元のEuclid空間にある点を保管するための空間分割データ構造である． $K$-d木を用いることによって，範囲探索や最近点探索などを高速に行うことができる．
$K$-d木は空間の情報を2分木として保管している． ここで2分木とは，根(root)とよばれる点（これは与えた点全体を表す）を2つの子点(child node) に分割する操作を繰り返すことによって生成されるデータ構造である．
2次元を例にすると，$K$-d木は以下のように構築される． 点集合を含む領域（長方形）を $x,y$軸のいずれか長い方の中央値で分割することによって2つの長方形を生成し，これらを子点とする． 子点には長方形に含まれる点を保持する． 各子点に対して同様の操作を繰り返し，長方形に含まれる点の数が一定数以下になったら，終了する． 子点をもたない点を葉(leaf)とよぶ．
コンストラクタの引数：
points: NumPyの形状 (npoints, ndim) をもつ浮動小数点数の配列； npointsは点数，ndimは次元数．
leafsize: 葉（子をもたない点）に含まれる点数を表す自然数；
メソッドquery(x)は，点（もしくは配列）xに含まれる各点からの最近点探索を行う． ここでx（の各要素）は，$K$-d木と同じ次元をもつ点の座標を表す配列とする． 返値は点xから近い順に並べ替えられた距離（の配列）と近い点の位置（を表す配列）である． オプションで与えることができる引数は以下の通り．
- k は近い点の数を表すパラメータであり，自然数を与える．既定値は $1$ である．
- p は距離を計算する際のノルムを表すパラメータである． $n$ 次元の点 $x,y$ 間のMinkowski $p$-ノルムは以下のように計算される． $$ \left( \sum_{i=1}^n |x_i -y_i|^p \right)^{1/p} $$
pは $1$ 以上の浮動小数点数を与え，$1$ のときは $L_1$（マンハッタン）ノルム，$2$ のときはEuclidノルム， 無限大のときは $L_{\infty}$ノルムになる．既定値は $2$ である．
例として， $[0,1]^2$ にランダムに発生させた $1000$個の2次元の点に対して， $(0.1,0.1)$ と $(0.5,0.5)$ から近い順に10個の点を表示する．
%matplotlib inline
from scipy.spatial import KDTree
import numpy as np
import matplotlib.pyplot as plt
points = np.random.rand(10000, 2)
tree = KDTree(points)
dis, near = tree.query(x=[[0.1, 0.1], [0.5, 0.5]], k=30)
print(dis) # 各点から近い点への距離（近い順）
plt.plot(points[:, 0], points[:, 1], "b+")
plt.plot(points[near, 0], points[near, 1], "ro")
None
メソッドquery_ball_point(x, r)は点（もしくは配列）xから距離r以内のすべての点のリスト（もしくはリストの配列）を返す．
その他の引数はqueryと同じである．
$[0,1]^2$ にランダムに発生させた $1000$個の2次元の点に対して， 中心 $(0.5,0.5)$ から $L_1$ノルムでの距離 $0.3$ 以内の点を表示を行う．
points = np.random.rand(1000, 2)
tree = KDTree(points)
ball = tree.query_ball_point(x=[0.5, 0.5], r=0.3, p=1)
plt.plot(points[:, 0], points[:, 1], "b+")
plt.plot(points[ball, 0], points[ball, 1], "ro")
None
凸包クラス ConvexHull
- 凸包：計算幾何学の基本となる概念；与えられたEuclid空間内の点を含む最小の凸多面体
コンストラクタの引数：
- points: NumPyの形状 (npoints, ndim) をもつ浮動小数点数の配列； npointsは点数，ndimは次元数．
凸包クラスConvexHullは以下の属性とメソッドをもつ．
- pointsは入力された点の座標を表す配列である．
- verticesは凸包を構成する点を表す配列である．2次元の場合には，反時計回りの順に保管されるが，3次元以上の場合には任意の順になる．
- simplicesは凸包を構成する面（単体）を保持する配列である． 面の数の配列であり，配列の各要素は点の次元の長さの配列である．
- neighborsは凸包を構成する面（単体）の隣接関係を表す配列である．配列の各要素は隣接する面の番号の配列である．
また，2次元の凸包を描画するための関数convex_hull_plot_2dも準備されている． これは凸包インスタンスを引数とし，matplotlibの図を返す．
from scipy.spatial import ConvexHull, convex_hull_plot_2d
points = np.random.rand(300, 2)
hull = ConvexHull(points)
convex_hull_plot_2d(hull)
print(hull.vertices)
Delaunay 3 角形分割クラス Delaunay
Delaunay3 角形分割: (2次元のときは）領域を点を通る3 角形に分割したとき，どの3角形の頂点を通る円も，他の点を含まない3角形分割
近い点同士を見つけるのに便利；Voronoi図と双対の関係
コンストラクタの引数：
- points: NumPyの形状 (npoints, ndim) をもつ浮動小数点数の配列； npointsは点数，ndimは次元数．
Delaunay3角形分割クラスDelaunayは以下の属性とメソッドをもつ．
- pointsは入力された点の座標を表す配列である．
- simplicesはDelaunay3角形分割に含まれる単体（2次元のときは3角形）を保持する配列である． 単体の数の配列であり，各要素は次元 $+1$ の長さの配列に単体の頂点番号が保管される．
- neighborsはDelaunay3角形分割に含まれる単体の隣接関係を表す配列である．配列の各要素は， 単体内の各点の反対側に隣接する単体の番号の配列である．対応する単体がない場合には $-1$ が保管される．
また，2次元のDelaunay3角形分割を描画するための関数delaunay_plot_2dも準備されている． これはDelaunay3角形分割インスタンスを引数とし，matplotlibの図を返す．
from scipy.spatial import Delaunay, delaunay_plot_2d
np.random.seed(5)
points = np.random.rand(80, 2)
tri = Delaunay(points, furthest_site=False)
fig = delaunay_plot_2d(tri)
Voronoi図クラス Voronoi
Voronoi図: 各点に近い空間で領域分けされた図
点に近い領域を見つけるのに便利；Delaunay 3 角形分割の双対
与えられた点を母点，各母点に近い空間から成る領域をVoronoi領域． Voronoi領域の境界をVoronoi境界，Voronoi境界の交点をVoronoi点
コンストラクタの引数：
- points: NumPyの形状 (npoints, ndim) をもつ浮動小数点数の配列； npointsは点数，ndimは次元数．
Voronoi図クラスVoronoiは以下の属性とメソッドをもつ．
- pointsは入力された点の座標を表す配列である．
- verticesはVoronoi点の座標を表す配列である．
- ridge_pointsはVoronoi境界を垂直に横切る母点対を表す． 実は，これはDelaunay3角形分割の枝（辺）に対応し，これがVoronoi図とDelaunay3角形分割が互いに双対といわれる所以である．
- ridge_verticesはVoronoi境界を表す点対である． Voronoi点の番号，もしくは点がVoronoi図の外側にある場合には $-1$ を保管する．
- regionsはVoronoi領域を表す点の配列である．点がVoronoi図の外側にあるときには $-1$ を保管する．
- point_regionは各母点が含まれるVoronoi領域の番号を表す．
また，2次元のVoronoi図を描画するための関数voronoi_plot_2dも準備されている． これはVoronoi図インスタンスを引数とし，matplotlibの図を返す．
例として格子上に配置した点に対するVoronoi図を求め，それを描画する．
同時に，対応するDelaunay3角形分割も示す．
from scipy.spatial import Voronoi, voronoi_plot_2d
points = np.array(
[[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]]
)
vor = Voronoi(points)
fig = voronoi_plot_2d(vor)
print("vertices=", vor.vertices)
print("reige_points", vor.ridge_points)
tri = Delaunay(points)
fig = delaunay_plot_2d(tri)
Delaunay3角形分割の例題に対するVoronoi図を計算し，対応するDelaunay3角形分割と比較する．
np.random.seed(5)
points = np.random.rand(8, 2)
vor = Voronoi(points)
fig = voronoi_plot_2d(vor)
tri = Delaunay(points) fig = delaunay_plot_2d(tri)
距離計算のためのサブモジュールdistance
euclidean(u,v): uとvの間のEuclid距離（直線距離）
$\|u-v\|_2 = \left(\sum_{i} (u_i-v_i)^2 \right)^{1/2}$
chebyshev(u,v): uとvの間のChebyshev距離
$(L_{\infty}$ ノルム) $\max_{i} |u_i -v_i|$
cityblock: uとvの間のManhattan距離（$L_{1}$ ノルム）
$\sum_{i} |u_i-v_i|$
minkowski(u,v,p): uとvの間のMinkowski距離（(Minkowski p-ノルム）
$\| u-v \|_p =\left(\sum_{i} |u_i-v_i|^p \right)^{1/p}$
from scipy.spatial.distance import euclidean
def NN_KD(points):
tree = KDTree(points)
near = {}
total_cost = 0.0
for i, p in enumerate(points):
dis, near_list = tree.query(p, k=2)
near[i] = near_list[1]
def NN_Voronoi(points):
vor = Voronoi(points)
dis, near = {}, {}
for i in range(len(points)):
dis[i] = np.inf
for (i, j) in vor.ridge_points:
d = euclidean(points[i], points[j])
if d < dis[i]:
dis[i] = d
near[i] = j
if d < dis[j]:
dis[j] = d
near[j] = i
n = 10000
points = np.random.rand(n, 2)
%time NN_KD(points)
%time NN_Voronoi(points)
import networkx as nx
def MST(points):
G = nx.Graph()
for i, p in enumerate(points):
for j, q in enumerate(points):
if i < j:
G.add_edge(i, j, weight=euclidean(p, q))
E = nx.minimum_spanning_edges(G)
def MST_Voronoi(points):
G = nx.Graph()
vor = Voronoi(points)
for (i, j) in vor.ridge_points:
G.add_edge(i, j, weight=euclidean(points[i], points[j]))
E = nx.minimum_spanning_edges(G)
return G, E
n = 200
points = np.random.rand(n, 2)
%time MST(points)
%time MST_Voronoi(points)
G, E = MST_Voronoi(points)
pos = {i: (points[i][0], points[i][1]) for i in range(n)}
nx.draw(G, pos=pos, node_size=10, edgelist=list(E))
lon = [
141.34694,
140.74,
141.1525,
140.87194,
140.1025,
140.36333,
140.46778,
140.44666999999998,
139.88361,
139.06083,
139.64889,
140.12333,
139.69167,
139.6425,
139.02361000000002,
137.21139,
136.62556,
136.22194,
138.56833,
138.18111000000002,
136.72222,
138.38306,
136.90667,
136.50861,
135.86833,
135.75556,
135.52,
135.18306,
135.83278,
135.1675,
134.23833,
133.05056000000002,
133.935,
132.45944,
131.47139,
134.55944,
134.04333,
132.76611,
133.53111,
130.41806,
130.29889,
129.87361,
130.74167,
131.6125,
131.42389,
130.55806,
127.68111,
]
lat = [
43.06417,
40.82444,
39.70361,
38.26889,
39.71861,
38.24056,
37.75,
36.34139,
36.56583,
36.39111,
35.85694,
35.60472,
35.68944000000001,
35.44778,
37.90222,
36.69528,
36.59444000000001,
36.06528,
35.66389,
36.65139,
35.39111,
34.97694,
35.18028,
34.73028,
35.00444,
35.021390000000004,
34.68639,
34.69139000000001,
34.68528,
34.22611,
35.50360999999999,
35.47222,
34.66167,
34.396390000000004,
34.18583,
34.06583,
34.34028,
33.84167,
33.55972,
33.606390000000005,
33.24944,
32.74472,
32.78972,
33.23806,
31.91111,
31.56028,
26.2125,
]
問題: Euclid最小費用完全マッチング （難）
2次元平面 $[0,1]^2$ 上にランダムに分布した偶数個の点に対する最小費用の完全マッチング （すべての点の次数が $1$ の部分グラフ）を求める問題を考える． ただし枝の費用は点間のEuclid距離とする．このような問題をEuclid最小費用マッチング問題とよぶ． $K$-d木を用いて，各点に対して近い順に $10$個の点の間にだけ枝をひいた疎なグラフを作成し， それ上で最小費用マッチングを求める近似解法を考える． $100$点のEuclid最小費用マッチング問題をランダムに生成し， 疎なグラフ上で求めた近似マッチングと， すべての点間の枝を用いて最適解を求めた場合の費用と計算時間を測定せよ．
ヒント： 最小費用マッチング問題の求解には，NetworkXの関数max_weight_matchingを用いることができる． この関数は負の重みには対応していない． したがって，費用の合計を最小化するためには，大きな数（たとえば枝の費用の最大値）から枝の費用を減じた値を新たに枝の費用と 定義して最大化を行う必要がある．
例として平均 $100$，標準偏差 $10$ の正規分布 $N(100,10^2)$ の統計量を分布に付随するメソッドstatで表示してみる． 引数のmoments='mvsk'は，平均(mean)，分散(variance)，歪度(skewness)，尖度(kurtosis)を意味する．
正規分布インスタンスを生成するnormは，位置パラメータlocで平均を，尺度パラメータscaleで標準偏差を与える．
密度関数 $f(x)$ をもつ連続な確率変数に対して， $X$ の期待値（平均）を $$ E[X]= \int_{-\infty}^{\infty} x f(x) dx $$ と定義する．
確率変数 $X$ と正数 $k$ に対する $k$次モーメント(moment)は，$m=E[X]$ としたとき， $E[ (X-m)^k ]$ で定義される．
分散(variance) $Var(X)$ は，2次モーメント $E[ (X-m)^2]$ と定義される． 分散は分布のばらつきを表す尺度である． $\sqrt{Var(X)}$ を標準偏差(standard deviation)とよび $\sigma$ と記す．
歪度(skewness)は $E[ (X-m)^3]/\sigma^3$，尖度(kurtosis)は $E[ (X-m)^4]/\sigma^4$ で定義され， それぞれ分布の非対称性ならびに尖りを表す尺度である．
正規分布の歪度は $0$ であり左右対称である．歪度が正の分布は右側の裾野が長くなり，逆に負の分布は左側の裾野が長くなる．
正規分布の尖度は $0$（$3$とする場合もあるが，SciPyでは$0$と定義されている）であり， 尖度が正の分布は，正規分布と比べて山の頂上付近で尖って，裾野が長い形状であり， 尖度が負の分布は，頂上付近が丸く裾野が短い形状をもつ．
from scipy.stats import norm
print(norm(loc=100, scale=10).stats(moments="mvsk"))
n = norm(loc=100, scale=100)
print(n.cdf(80))
from scipy.stats import gamma
g = gamma(1, scale=10) # a=1 (shape). 1/gamma=10 (scale)
print(g.stats(moments="mv")) # mean =a/gamma, variance =1/gamma^2
from scipy.stats import norm
n = norm(loc=100, scale=10)
print("array=", n.rvs(size=5))
print("matrix=", n.rvs(size=(2, 3)))
sf: 生存関数(survival function) （$1-F(x)$）
（離散，連続）確率変数 $X$ の生存関数は $$ P(X > x) $$ と定義される．
分布関数が与えた値以下の確率を返すのに対して，生存関数は与えた値より大きい確率を返す． （離散確率変数の場合には，「以上」でなく「より大きい」確率であることに注意する必要がある．）
例: 平均 100，標準偏差 10 の正規分布にしたがう需要をもつ商品に対して，123個仕入れたときの品切れ確率
from scipy.stats import norm print( norm.sf(123,loc=100,scale=10) )0.0107241100217
n = norm(loc=24, scale=8)
print(n.cdf(30))
print(n.ppf(0.9))
print(n.sf(30)) # もしくは 1- n.cdf(30)
print(n.ppf(1 - 0.1)) # 2と同じ
print(n.cdf(28) - n.cdf(20)) # 28以下の確率から20以下の確率を引く
print(24 - n.ppf(0.05)) # 確率が0.05（10%の半分）になるkを求め，平均24から引く
print(norm.expect(lambda x: abs(x)))
例題： 損出関数
在庫理論では，品切れ量を計算するために確率変数 $x$ が $k$以上になるときの期待値を 用いる．これは損出関数(loss function)とよばれ，以下のように定義される． $$ G(k)= \int_{k}^{\infty} (x-k) f(x) dx $$
たとえば，$k=0.1$ のときの標準正規分布の損出関数は， 期待値の計算する際の下限lbを $0.1$ としたときの $x-k$ の期待値として以下のように計算される．
損出関数は，密度関数 $f(x)$ と生存関数 $1-F(x)$ を用いて，解析的に $G(k)= f(k) -k(1-F(k))$ と計算できる．
以下のコードから， 結果が一致していることが確認できる．
k = 0.1
print(norm.expect(lambda x: x - k, lb=k))
print(norm.pdf(k) - k * norm.sf(k))
from scipy.stats import uniform
u = uniform(loc=0, scale=10)
print("平均と分散 ", u.stats(moments="mv"))
x = np.linspace(-3, 15, 100)
y = u.pdf(x)
plt.plot(x, y)
expon: 指数分布(exponential distribution)
指数分布とは，正のパラメータ $\lambda$ に対して，密度関数が $$ f(x) = \lambda e^{-\lambda x} \quad \forall x \geq 0 $$ で与えられる分布である．引数として指定する尺度パラメータscaleは $1/\lambda$ である． 分布関数は $$ F(x) = 1 - e^{-\lambda x} \quad \forall x \geq 0 $$ となる．平均は $1/\lambda$，分散は $1/\lambda^2$ である．
指数分布は無記憶性(memoryless property)をもつ唯一の連続分布 （離散分布では，後述する幾何分布が無記憶性をもつ．） として知られ， 稀な現象の時間間隔として用いられる
from scipy.stats import expon
ex = expon(scale=10)
print("平均と分散 ", ex.stats(moments="mv"))
x = np.linspace(-3, 100, 100)
y = ex.pdf(x)
plt.plot(x, y)
cauchy: Cauchy分布(Cauchy distribution)
Cauchy分布とは，分布の最頻値を表すパラメータ $x_0$， と広がりを表すパラメータ $\gamma$ に対して，密度関数が $$ \frac{1}{\pi\gamma\,\left[1 + \left(\frac{x-x_0}{\gamma}\right)^2\right]} $$ で与えられる分布である．
分布関数は， $$ \frac{1}{\pi} \arctan\left(\frac{x-x_0}{\gamma}\right)+\frac{1}{2} $$ である．Cauchy分布は平均をもたず，2次モーメントが無限大になる分布として知られている．
試しにSciPyで平均，分散，歪度，尖度を計算してみる．
from scipy.stats import cauchy
cau = cauchy(0, 1)
print(cau.stats(moments="mvsk"))
x = np.linspace(-10, 10, 100)
y = cau.pdf(x)
plt.plot(x, y)
平均，分散，歪度，尖度ともに，数値ではないことを表すnan(not a number)が表示された．
norm: 正規分布(normal distribution)
正規分布 $N(\mu,\sigma^2)$ は最も良く使われる連続分布であり， その密度関数は平均を表す位置パラメータ $\mu$ と標準偏差を表す尺度パラメータ $\sigma$ を用いて，
$$ f(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2} \right) $$と記述される．
中心極限定理から（標準偏差をもつ）独立な分布の平均の分布は正規分布に近づくことがいえる． 特に，平均 $0$，標準偏差 $1$ の正規分布 $N(0,1)$ は標準正規分布とよばれ，その密度関数と分布関数は，それぞれ $$ \phi (x) = \frac{1}{\sqrt{2\pi}} \exp\!\left(-\frac{x^2}{2} \right) $$ と $$ \Phi (x) = \int_{-\infty}^x \phi (t) dt = \frac{1}{2} +\frac{1}{2} erf \left( \frac{x}{\sqrt{2}} \right) $$ で定義される．ここで $erf(x)$ は誤差関数であり， $$ erf(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt $$ である．
lognorm: 対数正規分布(log-normal distribution)
対数をとったときに正規分布になる分布である． 正規分布にしたがう確率変数が負の値も許すのに対して，対数正規分布にしたがう確率変数は正の値しかとらないという性質をもつ． また，独立同分布に従う確率変数の積は漸近的に対数正規分布にしたがうという，大数の法則の積バージョンの性質をもつ．
密度関数は位置パラメータ $\mu$ と形状パラメータ $\sigma$（正規分布と異なり尺度パラメータでないことを注意）を用いて， $$ f(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}} \exp\!\left(-\frac{(\ln x-\mu)^2}{2\sigma^2} \right) \quad \forall x >0 $$ と記述される．分布関数は標準正規分布の分布関数を $\Phi$ としたとき，$\Phi \left(\frac{\ln x - \mu}{\sigma} \right)$ となる．
平均は $e^{\mu+\frac{ \sigma^2 }{2} }$，分散は $e^{2\mu+\sigma^2}(e^{\sigma^2}-1)$ である．
from scipy.stats import lognorm
ln = lognorm(1, loc=5)
print(ln.stats(moments="mv"))
x = np.linspace(-1, 10, 100)
y = ln.pdf(x)
plt.plot(x, y)
halfnorm: 半正規分布(half-normal distribution)
$X$ が正規分布 $N(0,\sigma^2)$ のとき， $Y=|X|$ がしたがう分布を半正規分布とよぶ． 一般に平均が $0$ でない場合には，折り畳み正規分布 (folded-normal distribution)とよばれるが， SciPyでは両方とも halfnorm で生成できる．
密度関数は位置パラメータ $\mu$ と尺度パラメータ $\sigma$ を用いて， $$ f(x)= \frac{1}{\sqrt{2\pi\sigma^2}} \, e^{ -\frac{(x-\mu)^2}{2\sigma^2} } + \frac{1}{\sqrt{2\pi\sigma^2}} \, e^{ -\frac{(x+\mu)^2}{2\sigma^2} } $$ と記述される．分布関数は標準正規分布の誤差関数を $erf(x)$ としたとき， $$ \frac{1}{2}\left[ erf\left(\frac{x+\mu}{\sqrt{2\sigma^2}}\right) + erf\left(\frac{x-\mu}{\sqrt{2\sigma^2}}\right)\right] $$ となる．
平均は $$ \mu_Y = \sigma \sqrt{\frac{2}{\pi}} \,\, \exp\left(\frac{-\mu^2}{2\sigma^2}\right) + \mu \, erf\left(\frac{\mu}{\sqrt{2\sigma^2}}\right) $$ であり， 分散は $\mu^2 + \sigma^2 - \mu_Y^2$ である．
from scipy.stats import halfnorm
hn = halfnorm(loc=0, scale=np.sqrt(2) * 10)
print(ln.stats(moments="mv"))
x = np.linspace(-1, 100, 100)
y = hn.pdf(x)
plt.plot(x, y)
erlang: Erlang分布(Erlang distribution)
Erlang分布は，独立で同一の $n$ 個の指数分布の和の分布であり， パラメータ $\lambda>0$ と自然数 $n$ に対して，密度関数が $$ f(x)={\lambda^{n} x^{n-1} e^{-\lambda x} \over (n-1)!} \quad \forall x \geq 0 $$ で与えられる分布である．ここで $n$ は形状パラメータ，$1/\lambda$ は尺度パラメータである． 分布関数は $$ F(x)= 1 - \sum_{k=0}^{n-1} {(\lambda x)^{k} \over k!} e^{-\lambda x} \quad \forall x \geq 0 $$ となる．平均は $n/\lambda$，分散は $n/\lambda^2$ である．
Erlang分布は待ち行列理論に応用をもつ．
from scipy.stats import erlang
er = erlang(3, scale=10)
print(er.stats(moments="mv"))
x = np.linspace(-1, 100, 100)
y = er.pdf(x)
plt.plot(x, y)
gamma: ガンマ分布(gamma distribution)
ガンマ分布は，正のパラメータ $\lambda, a$ に対して，密度関数が $$ f(x) = \frac{ \lambda^a x^{a-1} e^{-\lambda x}}{\Gamma(a)} \quad \forall x \geq 0 $$ で与えられる分布である．ここで $\Gamma$ はガンマ関数である． $a$ が形状パラメータ，$\lambda$ が尺度パラメータである． 分布関数は $$ F(x) = \frac{\gamma(a, \lambda x)}{\Gamma(a)} \quad \forall x \geq 0 $$ となる．ここで $\gamma$ は不完全ガンマ関数である． 平均は $a/\lambda$，分散は $a/\lambda^2$ である．
ガンマ分布はErlang分布のパラメータ $n$ を実数値 $a$ に一般化した分布である．
from scipy.stats import gamma
ga = gamma(3.5, scale=10)
print(ga.stats(moments="mv"))
x = np.linspace(-1, 100, 100)
y = ga.pdf(x)
plt.plot(x, y)
logistic: ロジスティック分布(logistic distribution)
ロジスティック分布は，パラメータ $\mu, s$ に対して，密度関数が $$ f(x) = \frac{\exp(-\frac{x-\mu}{s})}{s(1+\exp(-\frac{x-\mu}{s}))^2} $$ で与えられる分布である．分布関数は $$ F(x) = \frac{1}{1 + e^{-(x - \mu)/s}} = \frac{1}{2}\left\{\tanh\left(\frac{x - \mu}{2 s}\right) + 1 \right\} $$ となる．平均は $\mu$，分散は $\pi^s \mu^2/3$ である．
ロジスティック分布は正規分布と同様に釣鐘型の密度関数と$S$字（シグモイド）型の分布関数をもつが， ロジスティック分布の方が裾野が長いという特徴をもつ． ロジスティック分布の尖度は $1.2$ であり，正規分布の $0$ と比べて大きいので，より裾野が長いことが分かる．
from scipy.stats import logistic, norm
lo = logistic(0, 1)
print(lo.stats(moments="mvsk"))
x = np.linspace(-10, 10, 100)
y = lo.pdf(x)
y2 = norm(0, 1).pdf(x)
plt.plot(x, y)
plt.plot(x, y2)
weibull_min: Weibull分布(Weibull distribution)
Weibull分布は正のパラメータ $c, \lambda$ に対して，密度関数が $$ f(x) = c \lambda \left(\lambda x \right)^{c-1} \exp(-(\lambda x)^{c}) \quad \forall x \geq 0 $$ で与えられる分布である．ここで $c$ は形状パラメータ，$1/\lambda$ は尺度パラメータである． 分布関数は $$ F(x)=1-\exp( -(\lambda x)^c) $$ であり，平均は $\Gamma (1+1/c)/\lambda$，分散は $\left(\Gamma (1+2/c)-\Gamma (1+1/c)\right)^2/\lambda^2$ である． ここで $\Gamma$ はガンマ関数である．
Weibull分布は故障現象や寿命を記述するために用いられる． パラメータ $c<1$ のとき故障率が時間とともに小さくなり，$c>1$ のとき時間とともに大きくなる． $c=1$ のとき故障率は一定であり，この場合には指数分布と一致する．
from scipy.stats import weibull_min
wei = gamma(1, scale=10)
print(wei.stats(moments="mv"))
x = np.linspace(-1, 100, 100)
y = wei.pdf(x)
plt.plot(x, y)
beta: ベータ分布(beta distribution)
ベータ分布は正のパラメータ $a, b$ に対して，密度関数が $$ f(x)= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\, x^{a-1}(1-x)^{b-1} \quad \forall 0 \leq x \leq 1 $$ で与えられる分布である．ここで $\Gamma$ はガンマ関数であり，$a,b$ はともに形状パラメータである．
分布関数は $$ F(x)=I_x(a,b) $$ である．ここで $I_x$ は正則不完全ベータ関数である． 平均は $a/(a+b)$，分散は $ab/ \left( (a+b)^2 (a+b+1) \right)$ である．
ベータ分布は特定の区間の値をとる確率分布を表すときに有効であり，パラメータ $a,b$ を変えることによって様々な形状をとることができる． プロジェクト管理の手法として有名なPERT(Program Evaluation and Review Technique)においても確率的に変動する作業時間の近似として この分布を用いている．
from scipy.stats import beta
be = beta(2, 4)
print(be.stats(moments="mv"))
x = np.linspace(0, 2, 100)
y = be.pdf(x)
plt.plot(x, y)
t: $t$分布($t$ distribution)
$t$分布は正のパラメータ $\nu$ に対して，密度関数が $$ f(t) = \frac{\Gamma((\nu+1)/2)}{\sqrt{\nu\pi\,}\,\Gamma(\nu/2)} (1+t^2/\nu)^{-(\nu+1)/2} $$ で与えられる分布である．ここで $\Gamma$ はガンマ関数であり，$\nu$ は自由度(degree of freedom)とよばれる形状パラメータである． ちなみに $\nu$ の引数名はdegree of freedomの略でdfである．
分布関数は $$ F(t)=I_x(\nu/2,\nu/2) $$ である．ここで$I_x$ は正則不完全ベータ関数であり， $$ x = \frac{t+\sqrt{t^2+\nu}}{2\sqrt{t^2+\nu}} $$ とする． 平均は $0$，分散は $\nu>2$ のとき $\nu/(\nu-2)$，$1< \nu \le 2$ のとき無限大である．
いま $x_1,x_2,\ldots,x_n$ を平均 $\mu$ の独立な正規分布のサンプルとする． 標本平均と不偏分散を $$ \bar{x} = \frac{x_1+x_2+ \cdots+x_n}{n} $$ $$ s = \sqrt{ \frac{1}{n-1} \sum_{i=1}^n (x_i - \bar{x})^2 } $$ とする． このとき， $$ t = \frac{\bar{x} - \mu}{s/\sqrt{n}} $$ は，自由度 $\nu=n-1$ の $t$分布にしたがう． この性質を利用して，$t$分布は，標本値から母集団の平均値を統計的に推定する区間推定や，母集団の平均値の仮説検定に利用される．
平均との仮説検定はttest_1samp，独立な2つのサンプル間の仮説検定はttest_ind， 対応がある2つのサンプル間の仮説検定はttest_relで行うことができる．
ちなみに$t$分布は，$\nu=1$ のときCauchy分布になり，$\nu$ が無限大に近づくと正規分布に漸近する．
from scipy.stats import t
tdist = t(1)
print(lo.stats(moments="mvsk"))
x = np.linspace(-10, 10, 100)
y = tdist.pdf(x)
y2 = norm(0, 1).pdf(x)
plt.plot(x, y)
plt.plot(x, y2)
chi2: $\chi^2$分布($\chi^2$distribution)
$\chi^2$分布は正のパラメータ $k$ に対して，密度関数が $$ f(x)= \frac{1}{2^{\frac{k}{2}}\Gamma\left(\frac{k}{2}\right)}\; x^{\frac{k}{2}-1} e^{-\frac{x}{2}} \quad \forall x \geq 0 $$ となる分布である． ここで $\Gamma$ はガンマ関数であり， $k$ は自由度を表す形状パラメータである． $k$ の引数名はdfである．
分布関数は $$ F(x)= \frac{1}{\Gamma\left(\frac{k}{2}\right)}\;\gamma\left(\frac{k}{2},\,\frac{x}{2}\right) \quad \forall x \geq 0 $$ である．ここで $\gamma$ は不完全ガンマ関数である． 平均は $k$，分散は $2k$ である．
$x_1,x_2,\ldots,x_k$ を独立な標準正規分布のサンプルとするとき， $$ \sum_{i=1}^k (x_i)^2 $$ は自由度 $k$ の $\chi^2$分布にしたがう． この性質を利用して $\chi^2$分布は $\chi^2$検定chisquareやFriedman検定friedmanchisquareに用いられる．
from scipy.stats import chi2
ch = chi2(3)
print(ch.stats(moments="mv"))
x = np.linspace(0, 10, 100)
y = ch.pdf(x)
plt.plot(x, y)