Python言語による ビジネスアナリティクス 実務家のための最適化・統計解析・機械学習 科学技術計算モジュール SciPy

Presentation on theme: "Python言語による ビジネスアナリティクス 実務家のための最適化・統計解析・機械学習 科学技術計算モジュール SciPy"— Presentation transcript:

1 Python言語による ビジネスアナリティクス 実務家のための最適化・統計解析・機械学習 科学技術計算モジュール SciPy
久保幹雄

2 SciPyとは 様々な科学技術計算のための実装を含むモジュール群 bookpy/notebooks内で扱うモジュールは以下のとおり - optimize : 最適化と方程式の根（=解） - spatial : 計算幾何 - stats : 確率・統計 - interpolate : 補間 - integrate : 積分 - linalg : 線形代数 - cluster : クラスタリング - sparse : 疎行列

3 optimizeに含まれる 非線形最適化のための関数
関数 pythonでこの関数を実装する Google検索からのplot

4 optimizeに含まれる 非線形最適化のための関数
matplotlibとNumPyを用いてplotしてみる def f(x): 関数の定義 return x/(1+x**2) xx=np.linspace(-10,10,100) yy=f(xx) plot(xx,yy);None -10から10までを等間隔に 100等分した配列

5 optimizeに含まれる 非線形最適化のための関数
関数 の最小値を求める 関数 は1変数であるため minimize_scalarメソッドを用いる 引数methodでアルゴリズムの選択 - Brent : Brent法 - Golden : 黄金分割法 - Bounded : 最適解の範囲を引数boundsで与えた場合のBrent法

6 optimizeに含まれる 非線形最適化のための関数
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), method='Bounded')) アルゴリズムの選択 解の範囲（区間）

7 optimizeに含まれる 非線形最適化のための関数
問題1： 区間 における 関数 の最小値を求めよ fun: message: 'Solution found.' nfev: 27 status: 0 success: True x:

8 optimizeに含まれる 非線形最適化のための関数
複数の変数をもつ関数をある初期解からの探索で最小化 非凸な非線形関数の代表例：Rosenbrock関数

9 optimizeに含まれる 非線形最適化のための関数
minimizeで用いることができる解法　引数 : method Nelder-Mead法（非線形最適化用の単体法） - BFGS法（準ニュートン法） - 信頼領域Newton共役勾配法

10 optimizeに含まれる 非線形最適化のための関数
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') res = minimize(rosen, x0, jac=rosen_der, hess=rosen_hess, method='trust-ncg') 初期解（5次元） 解法によって必要な引数が 異なることに注意 勾配ベクトル ヘッセ行列

11 optimizeに含まれる 非線形最適化のための関数
問題2： Rosenbrock関数の最小化を初期値を変えて行え。 Rosenbrock関数は実行可能領域 で定義されている関数である。 この範囲外の初期値を与えるとどうなるか 実験せよ。

12 optimizeに含まれる 非線形最適化のための関数
制約付き最適化の例 半径　　　の円内　　　　　　　　　　の制約下で 線形関数　　　　　　を最大化 　　　　　　　　　線形関数　　　　　　　　　を最小化

13 optimizeに含まれる 非線形最適化のための関数
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)

14 optimizeに含まれる 非線形最適化のための関数
def f(x): x[0]をx , x[1]をyと考える 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) 　　　　　　　　

15 optimizeに含まれる 非線形最適化のための関数
def f(x): return -x[0]-x[1] bnds = ((0, None), (0, None)) f = lambda x: -x[0]-x[1] 無名関数lambdaで関数を定義 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)

16 optimizeに含まれる 非線形最適化のための関数
def f(x): return -x[0]-x[1] bnds = ((0, None), (0, None)) ((x下限, x上限),(y下限, x上 限)) 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)

17 optimizeに含まれる 非線形最適化のための関数
type : 関数が等式の場合は eq 非負の不等式の場合は ineq fun : 制約を表す関数を定義 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)

18 optimizeに含まれる 非線形最適化のための関数
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)

19 optimizeに含まれる 非線形最適化のための関数
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) fun : 最小化させたい関数 x0 : 初期値 (x,y) = (2,2) method : 解法 constraints : 制約

20 optimizeに含まれる 非線形最適化のための関数
問題3： 鳥が放物線 を描いて飛んでいる。 いま の位置にいるカメラマンが、鳥との距離が最小の 地点でシャッターを押そうとしている。鳥がどの座標に来たとき にシャッターを押せば良いだろうか？ 制約が等式であるため、 minimizeの引数methodはSLSQP SLSQP========== fun: jac: array([ , , ]) message: 'Optimization terminated successfully.' nfev: 28 nit: 7 njev: 7 status: 0 success: True x: array([ , ])

21 optimizeに含まれる 非線形最適化のための関数
大域的最適化と局所的最適化 - basinhopping(関数 , 初期解) 　焼きなまし法（アニーリング法） - brute(関数, 範囲） 　力づく法 大域的最適解 局所的最適解

22 optimizeに含まれる 非線形最適化のための関数
曲線へのあてはめ（最小2乗法） - 最小2乗法 : 測定で得られた数値の組の関係を近似した関数で 表すとき、残差の二乗和を最小とするような係数を決定し、近 似を行う方法 curve_fit(関数, data1, data2) * np.random.normal() 　標準正規分布に従う擬似乱数ベクトル

23 optimizeに含まれる 非線形最適化のための関数
曲線へのあてはめ（関数の根・不動点を求める） - brentq(関数, 値1, 値2) : Brent法によって関数の値1と値2の間に ある根を求める　* 関数の値の符号は異なる必要がある

24 Spatial 計算幾何学 - 幾何的問題 : 図形や空間に関する問題 spatialに含まれるクラス
計算機科学とは、幾何的問題を効率よく解くための基本アルゴリズム の体系化をめざす学問分野である　 - 幾何的問題 : 図形や空間に関する問題 spatialに含まれるクラス - K-d木 - 凸包 - Delaunay 3角形分割 - Voronoi図 - distanceモジュール

25 Spatial　計算幾何学 K-d木 K-d木（K-d Tree)とはK次元の平面空間にある点を保管するため の空間分割データ構造である。二分木の一種で根とよばれる点 を2つの子点に分割する操作を繰り返すことで得られる。

26 Spatial　計算幾何学 K-d木 K-d木（K-d Tree)とはK次元の平面空間にある点を保管するため の空間分割データ構造である。二分木の一種で根とよばれる点 を2つの子点に分割する操作を繰り返すことで得られる。 右図は空間を中央値で 分割するとき、 最初の操作後の平面空間を 二分木で表したもの

27 0 〜 1の一様乱数をもつ二次元ベクトルを1000個含む配列
Spatial　計算幾何学 points = np.random.rand(1000, 2) tree = KDTree(points) ball = tree.query_ball_point([0.5, 0.5], 0.1) plot(points[:,0], points[:,1], 'b+') plot(points[ball,0], points[ball,1], 'ro'); None 0 〜 1の一様乱数をもつ二次元ベクトルを1000個含む配列

28 Spatial　計算幾何学 points = np.random.rand(1000, 2) tree = KDTree(points) ball = tree.query_ball_point([0.5, 0.5], 0.1) plot(points[:,0], points[:,1], 'b+') plot(points[ball,0], points[ball,1], 'ro'); None KDTreeクラスのインスタンスオブジェクトは配列pointsを引数に要する 点 (0.5,0.5)から距離0.1以内の すべての点のインデックスの 配列を返す

29 Spatial　計算幾何学 points = np.random.rand(1000, 2) tree = KDTree(points) ball = tree.query_ball_point([0.5, 0.5], 0.1) plot(points[:,0], points[:,1], 'b+') plot(points[ball,0], points[ball,1], 'ro'); None 配列pointsに含まれるすべての要素の0番目

30 Spatial　計算幾何学 points = np.random.rand(1000, 2) tree = KDTree(points) ball = tree.query_ball_point([0.5, 0.5], 0.1) plot(points[:,0], points[:,1], 'b+') plot(points[ball,0], points[ball,1], 'ro'); None 配列ballに含まれるインデックスと同じインデックスもつ 配列pointsの要素の0番目

31 Spatial　計算幾何学 出力される画像

32 Spatial　計算幾何学 凸包（Convex hull）とは、 与えられたn次元Euclid空間内の点を含む最小の凸多面体 n次元の凸包はn-1次元の面から構成され、 構成された図形は単体（simplex）になる

33 Spatial　計算幾何学 points = np.random.rand(30, 2) hull = ConvexHull(points) convex_hull_plot_2d(hull) print( hull.vertices )

34 0 〜 1の一様乱数をもつ二次元ベクトルを30個含む配列
Spatial　計算幾何学 points = np.random.rand(30, 2) hull = ConvexHull(points) convex_hull_plot_2d(hull) print( hull.vertices ) 0 〜 1の一様乱数をもつ二次元ベクトルを30個含む配列

35 Spatial　計算幾何学 points = np.random.rand(30, 2) hull = ConvexHull(points) convex_hull_plot_2d(hull) print( hull.vertices ) ConvexHullクラスは コンストラクタの引数に n次元の点の配列pointsを要する

36 Spatial　計算幾何学 points = np.random.rand(30, 2) hull = ConvexHull(points) convex_hull_plot_2d(hull) print( hull.vertices )

37 Spatial　計算幾何学 points = np.random.rand(30, 2) hull = ConvexHull(points) convex_hull_plot_2d(hull) print( hull.vertices ) 端点

38 Spatial　計算幾何学 ボロノイ図（Voronoi diagram)とは、 各点に近い空間で領域わけされた図である。（双対=ドロネー図） 与えられた点を母点とよび、各母点に近い空間からなる領域を ボロノイ領域とよぶ。

39 Spatial　計算幾何学 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( vor.vertices )

40 Spatial　計算幾何学 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( vor.vertices ) 8個の点(x,y)を含む 配列pointsをnp.array で生成

41 生成した点データ（配列）を コンストラクタの引数に
Spatial　計算幾何学 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( vor.vertices ) 生成した点データ（配列）を コンストラクタの引数に

42 Spatial　計算幾何学 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( vor.vertices )

43 Spatial 計算幾何学 問題4： 右下図と同じボロノイ図をplotしなさい
from scipy.spatial import Voronoi, voronoi_plot_2d points = np.array([[1, 0], [1, 2], [2, 0], [2, 2]]) vor = Voronoi(points) fig = voronoi_plot_2d(vor) print( vor.vertices )

44 Spatial　計算幾何学 問題5： Weber問題（教科書88p 6.1 最適化より）で与えられている表6.1中の 家の位置をplotし、各家に近い領域で分けられたボロノイ図を描画せよ。 本問題　問1を解け。 本問題　問2を解け。 SLSQP========== [ ] [ ]

45 第1週はここまで

46 確率・統計 stats 確率・統計に関するモジュール stats 確率分布の基礎 共通のメソッド 連続確率変数 離散確率変数
データのあてはめ 相関と回帰 分布テスト

47 確率・統計 stats 正規分布クラス norm コンストラクタの引数 loc : 位置を表すパラメータ（μ : 平均） scale : 尺度を表すパラメータ（σ : 標準偏差） 既定値をもち、それぞれloc = 0, scale = 1 つまり標準正規分布 である

48 平均100,標準偏差10の正規分布 の確率分布関数 statsメソッドにより 統計量（平均,分散,歪度,尖度）を得る
from scipy.stats import norm print( norm(loc=100,scale=10 ).stats(moments='mvsk') ) 平均100,標準偏差10の正規分布 の確率分布関数 statsメソッドにより 統計量（平均,分散,歪度,尖度）を得る

49 確率・統計 stats 正規分布の固定化（変数に代入） n = norm(loc=100,scale=10) print( n.cdf(80) , n.pdf(80)) cdf (Cumulative density function : 累積分布関数)

50 確率・統計 stats 正規分布の固定化（変数に代入） n = norm(loc=100,scale=10) print( n.cdf(80), n.pdf(80) ) pdf（Probability density function : 確率密度関数）

51 確率・統計 stats 右図をプロットするプログラムコード from scipy.stats import norm import numpy as np n = norm(loc=100,scale=10) x = np.linspace(50,150,1000) y = n.pdf(x) plt.plot(x,y,label=r"$N(\mu=100,\sigma=10)$") plt.legend() plt.ylim((0,0.06)) plt.xlabel("x") plt.ylabel(r"$\varphi_{\mu,\sigma^2}(x)$")

52 確率・統計 stats 問題1： ある研究室でのコピー用紙の消費量は 需要200, 標準偏差120の正規分布にしたがうことがわかった。
コピー用紙の消費量の確率分布をプロットし確認しなさい コピー用紙の需要が80枚以下になる確率を求めなさい *プロットしたグラフをみて需要がどのような分布になっているか 考えて答えなさい

53 確率・統計 stats ガンマ分布 連続確率分布の一種で、2つの正のパラメータにより定義される。
形状パラメータ a：分布によって異なる。a=1のとき定義される分布は 指数分布となる。 尺度パラメータ scale

54 確率・統計 stats 確率変数オブジェクト（正規分布norm等）がもつメソッド rvs : 確率変数にしたがう擬似乱数を生成する - size : 乱数の数や配列の型を指定する n = norm(loc=100,scale=10) print('array=', n.rvs(size=5) ) print('matrix=', n.rvs(size=(2,3)) ) [ ]

55 確率・統計 stats 確率変数オブジェクト（正規分布norm等）がもつメソッド rvs : 確率変数にしたがう擬似乱数を生成する - size : 乱数の数や配列の型を指定する n = norm(loc=100,scale=10) print('array=', n.rvs(size=5) ) print('matrix=', n.rvs(size=(2,3) ) ) [[ ] [ ]]

56 右歪Gumbel分布： 右側の裾野が長い確率分布 確率変数の最大値が漸近的に従う分布
確率・統計 stats 確率変数オブジェクト（正規分布norm等）がもつメソッド stats : 平均、分散、歪度、尖度を返す。それぞれの頭文字 mvsk を引数にす ることで返り値を決定する。 g = gumbel_r(loc=0,scale=1) print( g.stats(moments='mvsk') ) x = np.linspace(-3, 5, 100) y = g.pdf(x) plot(x,y); None 右歪Gumbel分布： 右側の裾野が長い確率分布 確率変数の最大値が漸近的に従う分布

57 右歪Gumbel分布： 右側の裾野が長い確率分布 確率変数の最大値が漸近的に従う分布
確率・統計 stats 確率変数オブジェクト（正規分布norm等）がもつメソッド stats : 平均、分散、歪度、尖度を返す。それぞれの頭文字 mvsk を引数にす ることで返り値を決定する。 g = gumbel_r(loc=0,scale=1) print( g.stats(moments='mvsk') ) x = np.linspace(-3, 5, 100) y = g.pdf(x) plot(x,y); None 右歪Gumbel分布： 右側の裾野が長い確率分布 確率変数の最大値が漸近的に従う分布

58 確率・統計 stats 期待値：確率を乗じた平均 分散：分布のばらつきを表す尺度 標準偏差：分散の平方根
歪度：分布の非対称性を表す尺度 - 歪度が正の分布は右側の裾野が長く、負の分布は左側の裾野が長くなる 尖度：分布の尖りを表す尺度 - 尖度が正の分布は正規分布に比べて山の山頂付近で尖って、裾野が長く、　　　　 負の分布は頂上付近が丸く裾野が短い形状をもつ

59 平均100、標準偏差10の正規分布にしたがう 需要をもつ商品に対しての品切れ率を 5％にするための在庫量
確率・統計 stats 確率変数オブジェクト（正規分布norm等）がもつメソッド ppf : パーセント点関数（確率分布関数の逆関数） print( norm.ppf (0.95, loc=100, scale=10) ) 平均100、標準偏差10の正規分布にしたがう 需要をもつ商品に対しての品切れ率を 5％にするための在庫量

60 確率・統計 stats 確率変数オブジェクト（正規分布norm等）がもつメソッド sf : 生存関数（ ） print( norm.sf (123, loc=100, scale=10) ) 平均100、標準偏差10の正規分布にしたがう需要をもつ商品に対して、123個仕入れた時の品切れ確率を求める 1 - norm.cdf(123,loc=100,scale=10)と同値

61 確率・統計 stats 問題2： ある人気ラーメン屋では 平日の売り上げは平均150、標準偏差30の売り上げがあり
ある人気ラーメン屋では 平日の売り上げは平均150、標準偏差30の売り上げがあり 休日の売り上げは平均200、標準偏差20の売り上げがある。 平日の品切れ率を5%にするためには何杯分のスープの在庫を仕込んでお くべきか。 休日の品切れ率を3%にするためには何杯分のスープの在庫を仕込んでお くべきか。

62 120個以上の需要xが発生したとき、x-120個分の品切れ費用がかかる
確率・統計 stats 確率変数オブジェクト（正規分布norm等）がもつメソッド 平均100、標準偏差10の正規分布にしたがう需要をもつ商品に対して120個仕入れたとき、 品切れ費用が1個当り100円、在庫費用が1個当り10円としたときの期待費用を求める def f(x): if x>120: return 100*(x-120) elif x<120: return 10*(120-x) print( norm.expect(f, loc=100, scale=10) ) 120個以上の需要xが発生したとき、x-120個分の品切れ費用がかかる

63 120個以下の需要xが発生したとき、120-x個分の在庫費用がかかる
確率・統計 stats 確率変数オブジェクト（正規分布norm等）がもつメソッド 平均100、標準偏差10の正規分布にしたがう需要をもつ商品に対して120個仕入れたとき、 品切れ費用が1個当り100円、在庫費用が1個当り10円としたときの期待費用を求める def f(x): if x>120: return 100*(x-120) elif x<120: return 10*(120-x) print( norm.expect(f, loc=100, scale=10) ) 120個以下の需要xが発生したとき、120-x個分の在庫費用がかかる

64 確率・統計 stats 確率変数オブジェクト（正規分布norm等）がもつメソッド expect : 引数として与えられた関数funcの分布に対する期待値を求める。 - func : 関数 lb : 下限 ub : 上限 def f(x): if x>120: return 100*(x-120) elif x<120: return 10*(120-x) print( norm.expect(f, loc=100, scale=10) )

65 確率・統計 stats 問題3（発展問題）： ある人気ラーメン屋では 平日の売り上げは平均150杯、標準偏差30の売り上げがあり
ある人気ラーメン屋では 平日の売り上げは平均150杯、標準偏差30の売り上げがあり 休日の売り上げは平均200杯、標準偏差20の売り上げがある。 また、 ラーメンは1杯730円 スープの在庫は器材の関係により１５０杯分しか持つことができない 在庫を超える売り上げ数が見込まれたときは機会損失として考え、 そのときの機会損失額は１杯あたり1000円とする。 このラーメン屋の平日・休日の機会損失額を含めた売り上げ高の期待値を求めよ

66 確率・統計 stats SciPyでは正規分布クラスnorm以外にも様々な確率分布の確率変数オブジェ クトを使うことができるが、本スライドでは省略する。

67 確率・統計 stats 3次元プロットについてはNumPyを参考に
多変量正規分布クラス multivariate_normal ax = Axes3D(figure()) s = .1 X = np.arange(-3, 3.+s, s) Y = np.arange(-3, 3.+s, s) X, Y = np.meshgrid(X, Y) pos = np.dstack((X, Y)) n = multivariate_normal([0.0, 0.0], [[1.0, -0.8], [-0.8, 1.0]]) n = multivariate_normal([0.0, 0.0]) Z =n.pdf(pos) ax.plot_surface(X, Y, Z, rstride=1, cstride=1, norm=LogNorm(), cmap=cm.jet, linewidth=0.2); None 3次元プロットについてはNumPyを参考に

68 平均100、標準偏差10の正規分布にしたがう擬似乱数を10個生成
確率・統計 stats 最尤推定：最尤推定とは与えられたデータから確率分布のパラメータを推 定する方法である。 rvs = norm(loc=100,scale=10).rvs(size=10) print( norm.fit(rvs) ) 平均100、標準偏差10の正規分布にしたがう擬似乱数を10個生成

69 確率・統計 stats 最尤推定：最尤推定とは与えられたデータから確率分布のパラメータを推 定する方法である。 rvs = norm(loc=100,scale=10).rvs(size=10) print( norm.fit(rvs) ) >> ( , ) 与えられたデータの標本数を増やすと 推定されたパラメータがどのように変化するか確認してみよう （推定された平均,推定された標準偏差）

70 返り値のタプルは（標準偏差,平均,決定係数）
確率・統計 stats probplot : データの分布を正規分布と比較する正規確率プロットを作成する。 データが正規分布にしたがっているかどうかを判断する基準線や決定係数 （相関係数）がプロットされる rvs = norm(loc=100,scale=10).rvs(size=10000) xy_pair, corr = probplot(rvs, plot=plt) print( corr ) 返り値のタプルは（標準偏差,平均,決定係数）

71 補間 interpolate 補間（interpolation）とは、特定の領域内に与えられた離散的なデータ点を もとに、同じ領域内の他の点を近似する連続関数を構築するための方法で ある。補間は内挿ともよばれ、与えられた領域外の点の近似を行うことを 外挿・補外とも呼ぶ。

72 補間 interpolate x = np.arange(-4,5) y = np.sin(x)/(1+x**2) plot(x, y, 'o'); None

73 補間 interpolate lagrange : ラグランジュ多項式による補間 （すべての点を通る多項式を一つ決定する） f = interpolate.lagrange(x, y) print(f) xnew = np.arange(-4,4, 0.01) ynew = f(xnew) plot(x, y, 'o', xnew, ynew, '-') ; None

74 補間 interpolate interpld : 1次元データに対する補間。引数kindで補間の種類を指定する f = interpolate.interpld(x, y, kind= 'cubic') xnew = np.arange(-4,4, 0.01) ynew = f(xnew) plot(x, y, 'o', xnew, ynew, '-'); None

75 補間 interpolate 補間の種類 線形補間：区分ごとに異なる線分でつないだ区分的線形関数 引数kindはlinear

76 補間 interpolate 補間の種類 最近点補間：最も近いデータ点の値をとる関数 引数kindはnearest

77 補間 interpolate 補間の種類 0次補間：データ　　　　　　　　　の値をとる関数 引数kindはzero

78 補間 interpolate 補間の種類 スプライン線形補間：与えられた複数の点を通る滑らかな曲線 引数kindはslinear

79 補間 interpolate 補間の種類 スプライン2次補間：与えられた複数の点を通る滑らかな曲線 引数kindはquadratic

80 補間 interpolate 補間の種類 スプライン3次補間：与えられた複数の点を通る滑らかな曲線 引数kindはcubic

81 補間 interpolate 例題： あなたはテーマパークからアトラクションの待ち時間の長さを依頼された。 待ち時間を30分おきに測定した結果は開園時刻を0としたPythonのリストと して、与えられている。開園から閉園までの任意の時刻における待ち時間 を推定するにはどのようにしたら良いだろうか。 *待ち時間は最後尾に並んでいる人がアトラクションの会場に入場するまで の時間と定義する

82 補間 interpolate y = np.array([80,90,90,100,110,110,90,90,90,90,90,80, 80,60,60,60,60,70,70,70,60,35,40,50,45]) x=np.arange(0,len(y)) f = interpolate.interpld(x, y, kind= 'cubic') xnew = np.linspace(0,len(y)-1,100) ynew = f(xnew) plot(x, y, 'o', xnew, ynew, '-', lw=1,ms=7); None データの入力

83 補間 interpolate 問題4（教科書129p問題12）： ある新製品の毎日の需要量を調べたところ（教科書内の）表のよ うになっていることがわかった。ただし表内の – は需要を調べる のを忘れたことを表す。適当な補間を用いて、情報のない日の需 要量を予測せよ

84 補間 interpolate 問題4（教科書129p問題13）： あるアトラクションの待ち時間の長さを行った。待ち時間を30分 おきに測定した結果は、開園時刻を0としたPythonのリストとし て以下のように与えられている。開園から閉園までの任意の時刻 における待ち時間をLagrange補間と3字スプライン補間を用いて 推定せよ。ただし、リスト内のNaNはアトラクションが休止中で あったことを表す。 *NaNをどのように扱うか考えよう

85 積分 integrate quad : 1重定積分を求める。引数は関数f, 積分区間[a,b] (積分における部分区間の上限limitはoptional) 返り値は積分値と誤差のタプル (total,error) = integrate.quad(f, 0,len(y)-1,limit=1000) print( total/len(y) ) 平均

86 積分 integrate f = lambda x: math.pi*x**2
底面の半径rが10cm、高さhが10cmの円錐の体積を求める f = lambda x: math.pi*x**2 print( integrate.quad(f, 0, 10) )

87 積分 integrate dblquad : 2重定積分を求める。積分範囲の下限と上限を引数に要する 中心までの距離 f=lambda x,y: math.sqrt( (x-0.5)**2 +(y-0.5)**2 ) g=lambda x: 0.0 h=lambda x: 1.0 print(integrate.dblquad(f,0,1,g,h)) print((math.sqrt(2)+math.log(1+math.sqrt(2)))/6.0) 数値でなく関数で与えなければいけないことに注意

88 理論値は「はじめての確率論」（近代科学社）を参考に
積分 integrate dblquad : 2重定積分を求める。積分範囲の下限と上限を引数に要する 中心までの距離 f=lambda x,y: math.sqrt( (x-0.5)**2 +(y-0.5)**2 ) g=lambda x: 0.0 h=lambda x: 1.0 print(integrate.dblquad(f,0,1,g,h)) print((math.sqrt(2)+math.log(1+math.sqrt(2)))/6.0) 理論値は「はじめての確率論」（近代科学社）を参考に

89 積分 integrate 問題5： 楕円 の面積 は以下の式で求められる 楕円 の面積を求めよ

90 積分 integrate 問題6 東京海洋大学 編入学試験より

91 積分 integrate 問題（発展問題）： 面積1の円状の都市に一様に住んでいる人たちが、円の中心にあ るスーパーに行く時の距離の期待値を求めよ。 ちなみに理論値は である

92 線形代数　linalg 逆行列：A をn 次正方行列とするとA の逆行列は、 AB = BA = I を満たすn次正方行列Bである A = np.array([[1,-1,-1],[-1,1,-1],[-1,-1,1]]) print(linalg.inv(A)) print( np.dot(A,linalg.inv(A) )) 行列と逆行列のドット積（内積） は単位行列

93 LinAlgError : singular matrix
正則行列ではない行列の逆行列を求めようとすると... B = np.array([[2,-1,-1],[-1,2,-1],[-1,-1,2]]) print(linalg.inv(B)) print(np.dot(B,linalg.inv(B))) LinAlgError : singular matrix

94 線形代数　linalg 問題7 a = 3 のとき以下の問いに答えよ。 H26 東京海洋大学 個別学力検査より

95 線形代数　linalg 線形方程式を解く 線形方程式を行列とベクトルを用いて表現し、linalg.solveで解く A = np.array([[1,3,-2],[2,2,-1],[3,-1,1]]) b = np.array([-2,1,5]) print(linalg.solve(A,b))

96 線形代数　linalg 問題8 次の連立方程式を解きなさい。

97 線形代数　linalg 行列式（determinant）：正方行列Aの行列式 det（A）を求める *行列式=0となる行列の実行例 A = np.array([[1,2,3],[4,5,6],[7,8,9]]) b = np.array([-2,1,5]) print(linalg.inv(A)) # エラー print(linalg.solve(A,b)) # エラー print(linalg.det(A))

98 線形代数　linalg 問題9 東京海洋大学 編入学試験より

99 線形代数　linalg ノルム：平面や空間上のベクトルの長さのこと。linalg.norm関数 によって様々なノルムを求めることができる 引数を変更した場合の、実行結果は教科書p137を参考に

100 線形代数 linalg 固有値と固有ベクトルを求める関数 linalg.eig : 固有値と固有ベクトルを求める
linalg.eigvals : 固有値のみを求める A = np.array([[1,1,1],[1,2,3],[1,4,9]]) (v,X) = linalg.eig(A) print(v) print(X) V = np.diag(v) B = np.dot(X,np.dot(V,linalg.inv(X))) print(B) v : 固有値 X : 固有ベクトル

101 線形代数　linalg 問題10 東京海洋大学 編入学試験より