Python による科学技術計算の概要
神嶌 敏弘

www.kamishima.net
2022-04-27 更新

1

このチュートリアルについて
プログラミング言語 Python

Perl，Ruby，PHP と同様のスクリプト言語

1991年2月に Guido van Rossum が0.9.0を発表

数値関連パッケージが充実しており，データ分析でよく使われている
チュートリアルの前提と目的「 Pythonの文法」と「統計・機械学習」の知識を前提

ライブラリを用いたデータ分析用アルゴリズム実装の概要
チュートリアルのトピック

主要パッケージと実行環境

機械学習アルゴリズムの実装
関連パッケージと情報源

2

データ分析分野での Python
Kdnuggets による利用ソフトウェアの投票結果 2019年版

https://www.kdnuggets.com/2019/05/poll-top-data-science-machine-learning-platforms.html[参考] Top 10 Data Science Tools Over Time : https://youtu.be/pKPaHH7hnv8

Python はデータ分析ソフトとして重要な位置を占めている3

データ分析における Python
長所他のデータ分析ツールと比較したときの長所と短所

インタープリタだが主要な計算はネイティブで高速に実行

普通のプログラミング言語
データ分析以外のAPIとの連携が容易，データ構造やメモリ管理の自由度が大きい，デバッガ・テストなどの機能が利用できる数値計算系のライブラリや，他の数値計算ソフトへのライブラリが充実している短所

普通のプログラミング言語
言の記述などが必要になるクラスなどプログラミングに関わる宣ライブラリは充実しているが R には及ばない
4目次
第Ⅰ部：主要パッケージと実行環境

主要パッケージ：NumPy/SciPy と scikit-learn
実行環境：環境構築とクラウド

第Ⅱ部：機械学習アルゴリズムの実装

クラスの実装：scikit-learn形式のクラス と unittest など
数値計算tips：NumPyの基本と便利な機能
数値計算プログラミングの注意点

第Ⅲ部：関連パッケージと情報源

関連パッケージ：IPython，scikit-statmodels，最適化，ベイズ推定，数式処理，高速化など
情報源：Pythonによるデータ分析についての資料や勉強会

5

第Ⅰ部：主要パッケージと実行環境

主要パッケージ
NumPy/SciPy と scikit-learn

6利用されているPythonパッケージ
GitHub: Numpy and Scipy are the most popular packages for machine learning projects
2019−01−24

https://venturebeat.com/2019/01/24/github-numpy-and-scipy-are-the-most-popular-packages-for-machinelearning-projects/
7

NumPyNumPy：多次元配列を効率よく，簡単に扱うためのライブラリ
ホームページ： http://www.numpy.org/

内部的にはネイティブコードで実行されるので演算は高速
柔軟な要素の参照（fancy indexing）

データの一括処理（ユニバーサル関数，ブロードキャスティング）
数値・論理演算，線形代数，初等関数，集約演算，乱数生成など
フーリエ変換や関数あてはめなどのやや高度な演算も可能

歴史：Matrix Object (1994)，Numeric (1995)，Numarray(2001) などのプロジェクトを2005年に Travis Oliphant が統合して誕生
詳しい歴史： http://www.slideshare.net/shoheihido/sci-pyhistory
8

SciPySciPy：NumPy より高度な数値演算処理

ビルドに Fortran も用いているので，SciPy がインストールされていれば，NumPy のフーリエ変換なども高速になる
低次の演算：物理定数，疎行列，特殊関数，確率分布など

高次の計算：非線形最適化，補間，数値積分，計算幾何など

歴史：2001年に複数のライブラリを統合して誕生

NumPy / SciPy を読み込むときの一般的な省略名
import numpy as np
import scipy as sp
9scikit-learn
scikit-learn：非深層学習系の機械学習ライブラリの代表
ホームページ： http://scikit-learn.org/stable/

多くのアルゴリズムを統一された API で利用可能
歴史：2007に Google Summer of Code のプロジェクトとして始
まり，2010年以降 INRIA のメンバーが加わって発展した
主要な機能

教師あり学習：一般化線形モデル，SVM，アンサンブル学習など
教師なし学習：クラスタリング，異常検出，密度推定など

モデル選択と評価：交差確認，評価指標，超パラメータ探索など
データ変換：次元削減，標準化など

その他：ベンチマークデータ読み込み，テストデータ生成など
10scikit-learn によるロジスティック回帰

scikit-learn：線形回帰の例

In [1]: from __future__ import (print_function,
division,
absolute_import,
unicode_literals)

線形回帰の利用例

サンプルファイル：linear_regression.ipynb
パッケージの読み込み
In [2]: import numpy as np
from sklearn import linear_model

区間 [0, 10] と [0, 5] 上の一様分布に従う，2次元の独立変数のサンプルを100個生成
In [3]: X = np.array([np.random.uniform(0, 10, 100), np.random.uniform(0, 5, 100)]).T
回帰式 y = 1 x1 + 3 x2 + 10で従属変数を生成
In [4]: y = np.dot(np.array([[1, 3]]), X.T).ravel() + 10.0 + np.random.normal(
0, 0.1, 100)**線形回帰用のクラスを生成する；データに依存しないパラメータを指定**
ここでは切片項を使う指定を行う
In [5]: clr = linear_model.LinearRegression(fit_intercept=True)**fitメソッドでデータをあてはめ**11ここでは切片項を使う指定を行う

scikit-learn：線形回帰の例

In [5]: clr = linear_model.LinearRegression(fit_intercept=True)**fitメソッドでデータをあてはめ**
In [6]: clr.fit(X, y)Out[6]: LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normaliz
e=False)
従属変数の係数を調べると，ほぼ元の回帰式の 1 と 3 になっている
In [7]: clr.coef_Out[7]: array([ 0.99932678,

2.9981098 ])**predictメソッドで推定したモデルに基づく予測**
入力 x = (1, 2) に対する予測値は，ほぼ元の回帰式から得られる値 17 になっている
In [8]: clr.predict([[1, 2]])
Out[8]: array([ 16.99000072])

https://scikit-learn.org/stable/auto_examples/index.html
このExamplesの他に，User Guideにも多数のサンプル
12

第Ⅰ部：主要パッケージと実行環境

実行環境
環境構築とクラウド13

実行環境の構築
Python バージョン3を使うこと

Pythonバージョン2のサポートは2020年1月1日でサポート終了
2019年以降リリースのNumPyはバージョン3のみサポート

https://github.com/numpy/numpy/blob/master/doc/neps/dropping-python2.7-proposal.rst

Python 2 の資産がある場合：six などの移行ライブラリなどを使ったり，2to3などのツールを使ってPython3移行
Python 3：統計量パッケージの標準化や行列演算子（3.5以降）やデータクラス（3.7以降）など数値計算に有用な機能

Python による科学技術計算の実行環境構築

pip コマンドを使ってパッケージをインストールするのは問題も多かったが，最近は改善されてきた各種のパッケージをまとめたインストーラは容易に使えるが，pipでインストールしたものと整合性を損なう場合もある

14

インストーラ
フリーのものと商用のものがあるが，商用でもパッケージ数やサポートに制限のある無料版が提供されている
PyMC などインストーラによっては含まれないものもフリー

apt / yum / port / brew（Linux / Mac / BSD などUNIX系）

商用
Continuum Analytics Anaconda
https://store.continuum.io/cshop/anaconda/
Enthought Canopy
https://www.enthought.com/
その他 http://ibisforest.org/index.php?python%2Fnumpy を参照
15

第Ⅱ部：機械学習アルゴリズムの実装

クラスの実装scikit-learn形式のクラス と unittest

16

アルゴリズムのクラスによる実装
機械学習のアルゴリズムは，関数でも実装できるが，クラスを定義して実装すると以下の利点がある
Pythonの適性：Python はオブジェクト指向型言語なので，クラスによる実装に適している
クラスの継承 や Mixin の利用：一部だけが異なる学習アルゴリズム
を容易に実装できる
学習結果の保存：cPickle といったオブジェクトのシリアライズを利用して，学習したモデルのオブジェクトを保存可能
scikit-learnとの連携：scikit-learn のクラスの作成規則に沿って実装することで，その機能を利用できる
17

scikit-learnのAPI仕様
scikit-learn のクラス設計の基本仕様
Coding Guidelines：http://scikit-learn.org/stable/developers/
コンストラクタ：クラスの初期化

引数：データに依存しないアルゴリズムのパラメータﬁt()メソッド：あてはめ・学習

引数：訓練データと，データに依存したパラメータ

predict()メソッド：学習済みデータを用いた予測
引数：予測対象の新規の入力データ

score()メソッド：モデルのデータへのあてはめの良さの評価
評価対象のデータを，このメソッドの引数で指定する．

transform()メソッド：次元削減などのデータ変換

18

MixIn を用いたクラス
Mixin（メソッドのみのクラス）による一部だけが異なる手法の実装
LinearFuncMixin
objectfuncメソッドは線

FuncFitMixin
objectfuncメソッドで，あてはめを形関数

するﬁtと予測をするpredict

LogFitBaseEstimator, RegressorMixin,
LogFuncMixin, FuncFitMixin
対数関数でのあてはめ

LinearFit
BaseEstimator, RegressorMixin,
LinerFuncMixin, FuncFitMixin
線形関数でのあてはめ

LogFuncMixin
objectfuncメソッドは対数関数

緑色は親クラス
BaseEstimator と RegressorMixin は scikit-learnのクラス
必要な機能を組み合わせて部分的に異なるクラスを容易に実装可能
19

ook

from scipy.optimize import curve_fit
from sklearn.base import BaseEstimator, RegressorMixin
import matplotlib.pyplot as plt

Mixin を用いたクラス

2014/11

In [1]: import numpy as np
In [2]: class LogFuncMixin(object):from scipy.optimize import curve_fit
def func(self, x, a, b):from sklearn.base import BaseEstimator, RegressorMixin
return a * np.log(x)+ bimport matplotlib.pyplot
as pltサンプルファイル：class̲mixin.ipynb
In [1]: import numpy as np
対数関数のMixinfrom scipy.optimize
import curve_fit
Inclass
LinearFuncMixin(object):
In [3]:[2]: class
LogFuncMixin(object):
fromdefsklearn.base
BaseEstimator,
RegressorMixindeffunc(self,importx, a,a, b):b):func(self,x,import matplotlib.pyplot
return aa ** np.log(x)x + b as plt
return+ b

Inclass LinearFuncMixin(object):
LogFuncMixin(object):
FuncFitMixin(object):[4]:
In [2]:[3]: classfunc(self,x,x,y):a, b):b):
def func(self,fit(self,defx,a,returnpcov
a ** x
np.log(x)
+ b
popt,
=
curve_fit(self.func,
x, y)returna+bself.coef_ = poptIn
class FuncFitMixin(object):LinearFuncMixin(object):
In [3]:[4]: classdef predict(self, x):
def fit(self,func(self,x,x,y):a, b):
def
return self.func(x,*self.coef_)
returnpcov
a * =x curve_fit(self.func,
+ b
popt,
x, y)
self.coef_ = popt
LinearFit(BaseEstimator, RegressorMixin,
In [5]: class
In [4]: class FuncFitMixin(object):LinearFuncMixin,
FuncFitMixin):
defx):
def predict(self,fit(self, x, y):passreturn self.func(x, *self.coef_)
popt, pcov = curve_fit(self.func, x, y)

self.coef_ = popt
RegressorMixin,
In [6]:
LogFit(BaseEstimator, RegressorMixin,[5]: class LinearFit(BaseEstimator,
def predict(self,x): FuncFitMixin):
LinearFuncMixin,
FuncFitMixin):
LogFuncMixin,
passreturn self.func(x, *self.coef_)

線形関数のMixin
あてはめと予測のMixin

対数関数へのあてはめ
線形関数へのあてはめ

LinearFit(BaseEstimator,
RegressorMixin,[6]: class
LogFit(BaseEstimator,
RegressorMixin,[5]:
In [7]:
x = np.linspace(0.1,
2.0, 100)
LinearFuncMixin,
FuncFitMixin):0.2, 100)
LogFuncMixin,
y = 2.0 * np.log(x)+ 10 + FuncFitMixin):
np.random.normal(0,pass

In [8]: clf1 = LogFit()[7]: x= np.linspace(0.1,
2.0, 100)
In [6]:
class
LogFit(BaseEstimator,
RegressorMixin,
clf1.fit(x,
y)
y
= 2.0
* np.log(x)+ 10 + np.random.normal(0,0.2, 100)
LogFuncMixin,
FuncFitMixin):
clf2= LinearFit()

20

scikit-learn を用いて，交差確認を行う

scikit-learnとの連携

http://scikit-learn.org/stable/modules/cross_validation.html (http://scikitlearn.org/stable/modules/cross_validation.html) の例題
In [3]: from __future__ import (print_function,
division,
absolute_import,
unicode_literals)

scikit-learn のクラスの規則に従うことで，その機能を利用する
fit() や predict() などのメソッドの規則を踏襲する
sklearn.base.BaseEstimator クラスを親とし，score() メソッドをパッケージの import

含む sklearn.base.ClassifierMixin なども継承する
In [4]: import numpy as np
from sklearn import model_selection
from sklearn import datasetsfrom sklearn import svm

交差確認による評価や，超パラメータ最適化が容易になる
サンプルファイル：cross_validation.ipynb

IRISデータのロードと SVM 分類器の作成
In [5]: iris = datasets.load_iris()
clf = svm.SVC(kernel='linear', C=1)
交差確認による汎化誤差の評価． n_jobs=-1 とすると全CPUを使って並列計算
In [6]: scores = model_selection.cross_val_score(clf, iris.data, iris.target,
cv=5, n_jobs=5)
np.mean(scores)
Out[6]: 0.9800000000000001

21自動テストを書こう
unittest：プログラムの最小単位ごとに，その入出力を検証する

厳密に最小単位でなくても，リファクタリングや効率化で書き換えたときに自動で検証できるようになる
一度信頼できる結果を得たときに，その結果と一致するかどうかだ
けの簡単なものでも有用なので，テストを書くことを薦める

Python の場合は nosetests というコマンドによりテストを実行できる

パッケージ全体について，全てのテストをまとめて実行できる




22

unittestには何を書いたらいいのか？
本来は，エラー処理を含めて全てのコードが実行されるように書く
あくまで理想，変更が頻繁な実験コードなどでは無理がある
外部APIの仕様確認

機械学習のアルゴリズムを実装し，クラスを作って，ﬁt() メソッド
で学習し，predict() メソッドで予測する場合
小さなテストデータを読み込んで，学習後の損失関数値やパラメー
タ値，また，予測結果が一致するかを検証
デバッグのときに書いたコードに，実行結果を加えて簡単に作っておく

特徴や目的関数の値が全て同じなどの特殊な境界条件のデータも追加で検証しておくと安全23

外部APIの動作検証
サンプルファイル：test_svc.py


㵸‫׼־ذ٭ظعتطם׈‬‫⮆׾׌عتط‬걉㉺‫ظٓך‬ٜ؅㰢綗





➟丗ֿ┉舙‫־׾ַי׊‬؅吾錞㲔丗⡁‫⸬ע‬㳡‫┉ם‬舙‫ׂםעך‬劼‫־‬劼‫סַ׼׃‬瞿䈱‫┉ך‬舙؅吾錞


⮆걉磵冽ֿ┉舙‫־׾ַי׊‬؅吾錞‫ؠ‬ٚ‫םֹ׻סت‬ꦕ丆⡁‫⸬ע‬㳡‫┉ם‬舙؅吾錞


24

unittestと乱数
乱数で初期化するアルゴリズムでは，実行するたびに結果が変わる
ため結果が一致するかどうかを検証できない
乱数のシードを与え，クラス専用の疑似乱数生成器を作成する
scikit-learn のユーティリティ関数 check_random_state() が便利
サンプルファイル：random_state.py

‫⮣סغ٭ب‬僿⴫؅车‫ַם؂‬ֿ‫ؚنظס‬ٜ‫ع‬⡁‫⫋׷‬杯䓪ֿ䑒釐‫ם‬㲔꽦‫עך‬鸵䎎‫ם‬丝丗⡁؅┙ֻ‫׾‬


ْ‫▞עךغشخ‬丗榟䡗㉺؅暠䐂





25

unittestには何を書いたらいいのか？
内部メソッドの検証

最適化する目的関数（損失関数）やその勾配関数などの値を検証

単純な別コードや手計算で求めた結果と比較できる数10個の非常に小さなデータを使って，計算結果を信頼できるように確率なら総和が 1 になるなど，結果が当然満たすべき条件も有用

パラメータの値が全て 0 や全て 1，関数値が 0 になる点など，境界条件や特異値を検証コードに入れておくパラメータや潜在変数の初期値は，やはり境界条件や一様といった特別な条件で試しておく

26

境界条件の例
pLSA_initialization

pLSA アルゴリズムは，EMアルゴリズムという交互最適化法の一種
In [6]: n = 100
k = 2
d = gen_data(n)を用いて解を求める# print(d)2018/04/30 21)18

pXgZ = np.tile(1 / n, (n, k))pYgZ = np.tile(1 / n, (n, k))pZ = np.tile(1 / k, k)
for i in range(5):pXgZ, pYgZ, pZ = pLSAstep(d, pXgZ, pYgZ, pZ)
print(log_loss(d, pXgZ, pYgZ, pZ))

教科書などでは初期値について，あまり記述はないが一様な値で初期化すると動かない
データが偶然に一様だったりすると動かないといった場合もパラメータを一様に初期化
一様な定数

pLSA_initialization

In [6]: n = 100
k = 2
d = gen_data(n)# print(d)pXgZ = np.tile(1 / n, (n, k))pYgZ = np.tile(1 / n, (n, k))pZ = np.tile(1 / k, k)
for i in range(5):pXgZ, pYgZ, pZ = pLSAstep(d, pXgZ, pYgZ, pZ)
print(log_loss(d, pXgZ, pYgZ, pZ))
-44014.18824-44014.18824-44014.18824-44014.18824-44014.18824

対数尤度は不変-44014.18824-44014.18824-44014.18824-44014.18824-44014.18824

パラメータを乱数で初期化

Dirichlet分布に従う乱数

パラメータ Pr[X|Z], Pr[Y|Z] , Pr[Z] を乱数（Dirichlet分布）で初期化した場合
2018/04/30 21)18

対数尤度は順調に増加

In [7]: n = 10k = 2
d = gen_data(n)# print(d)pXgZ = np.random.dirichlet(alpha=np.repeat(1 / n, n), size=k).T
pYgZ = np.random.dirichlet(alpha=np.repeat(1 / n, n), size=k).T
pZ = np.random.dirichlet(alpha=np.repeat(1 / k, k))for i in range(5):pXgZ, pYgZ, pZ = pLSAstep(d, pXgZ, pYgZ, pZ)
print(log_loss(d, pXgZ, pYgZ, pZ))
-209.233267896-202.685252472-198.022879607-192.019675469-186.921832146

対数尤度は増加

サンプルファイル：pLSA_initialization.ipynb

パラメータ Pr[X|Z], Pr[Y|Z] , Pr[Z] を乱数（Dirichlet分布）で初期化した場合
対数尤度は順調に増加

27

第Ⅱ部：機械学習アルゴリズムの実装

数値計算tips
NumPyの基本

28

NumPy 配列と多重リスト配列の違い
np.ndarray：N-d Array すなわち，N次元配列を扱うためのクラス

1. メモリ上での保持多重リスト：リンクでセルを結合した形式でメモリ上に保持
多重リストは動的に変更可能

np.ndarray：メモリの連続領域上に保持形状変更には全体の削除・再生成が必要．
2. 要素の型

多重リスト：リスト内でその要素の型が異なることが許される
np.ndarray：基本的に全て同じ型の要素

29

NumPy 配列と多重リスト配列の違い
3. 配列の形状

多重リスト：次元ごとに要素数が違ってもよい
np.ndarray：各次元ごとの要素数が等しい
行ごとに列数が異なるような2次元配列などは扱えない

4. 高度で高速な演算操作
np.ndarray：行や列を対象とした多くの高度な数学的操作を，多重リストより容易かつ高速に適用可能 (fancy indexing)

配列中の全要素，もしくは一部の要素に対してまとめて演算や関数を適用することで，高速な処理が可能 (ブロードキャスト，ユニバーサル関数)30

NumPy 配列の生成
np.array を使った生成
np.array(object, dtype=None)object には，配列の内容を，array̲like という型で与える
array̲like型：配列を np.ndarray の他，（多重）リストや（多重）タプルで表現したもの
要素が 1, 2, 3 である長さ 3 のベクトルの例




タプルを使った表現も可能です


31

NumPy 配列の生成
2重にネストしたリストで表した配列の例






リストの要素に np.ndarray タプルを含むことも可能





32

0行列と1行列
0行列（要素が全て 0 の配列）の生成
np.zeros(shape, dtype=None)
1行列（要素が全て 1 の配列）の生成
np.ones(shape, dtype=None)
shape：スカラーや，タプルによって配列の各次元の長さを表現
長さが 3 の 0 ベクトルの例



3 4の1行列の例（引数をタプルにすることを忘れないように）




33

初期化なし配列

要素の初期化なしで指定した大きさの配列を生成
np.empty(shape, dtype=None)配列の生成後，その内容をすぐ後で書き換える場合には，配列の要
素を 0 や 1 で初期化するのは無駄メモリだけ確保・初期化しな
い

34

単位行列
単位行列を生成np.identity(n, dtype=None)n は行列の大きさを表す

例：4 と指定＝大きさ4 4の行列（単位行列は正方行列なので）






35

NumPy 配列の属性と要素の参照
np.ndarray クラスの主な属性
class np.ndarraydtype：配列要素の型
後ほど，詳細を述べる
ndim：配列の次元数ベクトルでは 1 に，配列では 2
shape：配列の各次元の大きさ＝配列の形状
各次元ごとの配列の大きさをまとめたタプルで指定
例：長さが 5 のベクトルは (5,)
※ Python のタプルは1要素のときは “,” が必要
例：2 3行列では (2, 3)
36

NumPy の dtype
Python のビルトイン型に対応する以下の型はそのまま用いて大きな問題はない真理値型＝ bool，整数型＝ int，浮動小数点型＝ float，
複素数型＝ complex

これらの型に対する正式な NumPy の型は，np. を前に付けるだけ
np.bool, np.int, np.float, np.complex※ さらに最後に “̲” を最後につけた np.int_ などが，Pure Pythonと厳密には互換の型

メモリのビット数を明示的に表す型：np.int32 や np.float64 が定義されていて，次の目的で利用するメモリを特に節約したい場合（10までなら8bitで十分）

C や Fortran で書いた関数とリンクする（32bit・64bitのどちらの環境で実行しても暴走しない）
37

NumPy の dtype（文字列）
ビルトイン型の str / unicode と NumPy のdtype の相違点
np.ndarray の要素の大きさが同じである必要 → 文字列は固定長
NumPyでの文字列型：NumPy の型を返す関数 np.dtype() を利用
np.dtype('S<the length of string>')
例：最大長が16である文字列


Unicode文字列の場合は，この S が U に置き換わります

※ Python3 の文字列は Unicode型 として指定

38

配列のdtypeの指定方法
np.array() などの配列生成関数の dtype 引数で指定する方法







np.ndarray の np.ndarray.astype() メソッドを使う方法








39

np.ndarrayの要素の参照
多様な要素の参照方法がありますが，ここでは基本的な方法のみ
各次元ごとに何番目の要素を参照するかを指定する方法
例：1次元配列であるベクトル a の要素 3 を a[3] 参照




※ 添え字の範囲は，1 からではなく 0 から数える※ a.shape[0] で，第1次元の要素の長さとして 5 が得られたとき添え字の範囲はそれより一つ前の 4 まで





40

1次元配列のスライス
スライス：リストや文字列などのスライスと同様の方法により，配列の一部分をまとめて参照する方法
1次元配列：リストのスライス表記と同様の 開始:終了:増分 の形式※ 負値は最後からの個数を示す









41

複数要素の同時指定
配列やリストを使って複数の要素を指定し，それらをまとめた配列を作る
11行目の例：リストを使って，0番目，2番目，1番目，2番目の要
素を繋げた配列
12行目の例：配列を使って要素を取り出して，まとめた例








42

2次元配列のスライス
スライスは，2次元以上の配列でも同様の操作が可能
特に， “:” のみを使って，行や列全体を取り出す操作は頻繁に利用












43

np.ndarrayと数学の行列※ 1次元の np.ndarray 配列には，線形代数でいう縦ベクトルや横ベク
トルという区別はなく，1次元配列は転置できない

縦ベクトルや横ベクトルを区別して表現するには，それぞれ列数が1である2次元の配列と，行数が1である2次元配列を利用
縦ベクトルの例：







横ベクトルの例：



44

配列の次元数の追加
np.newaxis を利用して1次元のベクトルを2次元の行列に変換








reshape メソッドをしてベクトルを配列に変換※ -1 を指定すると，全体の要素数が不変となるように大きさを計算







45

第Ⅱ部：機械学習アルゴリズムの実装

数値計算tips
NumPy/SciPyの便利な機能

46ユニバーサル関数
ユニバーサル関数：入力した配列の各要素に関数を適用し，その結果を入力と同じ形の配列して返す
ユニバーサル関数の機能を利用するには，mathパッケージのmath.log() などではなく，NumPy の np.log() を用いる





※ 論理関数のユニバーサル関数は and などではなく np.logical_and()などを用いる47

ユニバーサル関数の作成
ユニバーサル関数ではないユーザ関数をユニバーサル化するには
np.vectorize() を用いる

np.vectorize(pyfunc)
pyfunc で指定したユーザ関数をユニバーサル化した関数を返す※ 関数を引数にとり，関数を出力するのでデコレータとして利用可能





※ 出力が複数ある関数をユニバーサル化する np.frompyfunc() も48ブロードキャスト
ブロードキャスト：要素ごとの演算を行うときに，配列の大きさが異なっていれば自動的に要素をコピーして大きさを揃える機能
例：大きさがそれぞれ 0次元配列 a，長さ3の1次元ベクトル b，
2 3 の行列 c の要素積を求める
0

a
1

2

0

0 00

00

00

1 00

00

00

b
1

2

0

0 00

01

02

1 00

01

02

c
1

2

0 00

01

02

1 10

11

12

大きさが1の次元は，必要に応じて要素の値を自動的にコピーする
np.newaxis により次元数を合わせてから利用することを薦める
49

行列演算の関数プログラミング的実装
map() と reduce() を使った関数プログラミング的な実装に基づくと，NumPy による計算は理解しやすいだろう
リストではなく，共通の配列に map() して要素ごとの計算をしたあと，sum() や max() などの集約演算で reduce() すると考える
例：pLSAモデル

Pr[X, Y ] = ∑Z Pr[X ∣ Z] Pr[Z ∣ Y ] Pr[Y ]Pr[X | Z ]pXgZ[x, z]

Pr[Z | Y ]pZgY[z, y]Pr[Y ]
pY[y]

z

y

y

x

z

50

行列演算の関数プログラミング的実装
np.newaxisやブロードキャストで共通の配列に揃える

pXgZ[x,z,np.newaxis]pZgY[np.newaxis,z,y] pY[np.newaxis,np.newaxis,y]zznewaxisx

newaxisynewaxisnewaxisy

要素ごとの計算（map のイメージ）

pXgZ[x,z,np.newaxis] * pZgY[np.newaxis,z,y] * pY[np.newaxis,np.newaxis,y]

集約演算（reduceのイメージ）

np.sum(pXgZ[x,z,np.newaxis] * pZgY[np.newaxis,z,y]* pY[np.newaxis,np.newaxis,y], axis=1)

z
x

sum()

y

y
x

51

for文のブロードキャストへの変換
for文を用いた計算を，ブロードキャストを使ったものと置き換える

 個のデータを含む配列  と  は，それぞれ  と の値をとる

 と  の対の値を数え上げ分割表を作る
for文による基本的な実装

 個のデータを走査し，一つずつ数え上げる

サンプルファイル：contingency_table_example.ipynb



※ 参考： http://www.kamishima.net/mlmpyja/nbayes2/distclass.html
52

for文のブロードキャストへの変換
配列を参照する変数が，for文の変数になるようにする
配列 ct のインデックス‫ס‬などを，for文の変数 xi に置き換え
 や  の値を番号ではなく 0/1 ベクトルの形式で表す目的 や  の値が，for文の変数  や  と一致したとき数え上げる
 の値のために新たなfor変数  を導入






53

for文のブロードキャストへの変換
ブロードキャストを用いた実装

各for変数を，配列の各次元に割り当てる

for変数 ，， をそれぞれ配列の次元0, 1, 2に割り当てる

for文の変動範囲を np.arange で作った配列と置き換える

さらに，3次元（for文の段数）の設定した次元に割り当てる

if 文の判定はユニバーサル関数でまとめて実行（map演算）
最後に次元0の変数 i について総和をとる（reduce演算）for変数  は次元0に割り当てた





54

欠損値や無限値
欠損値 np.nan や無限大・無限小 np.inf などを含むデータの処理
配列の要素が全て有限値かどうかを検査する

np.all(np.isfinite(x)) の代わりに np.isfinite(x.sum())

np.nan を除外した集約演算
総和 np.nansum() や最大値 np.nanargmax() などがある
この場合に特化した bottleneck というパッケージも存在
欠損値を他の値で置き換える

ユニバーサルな3項演算子にあたる np.where() が便利





55

欠損値や無限値
scikit-learn で開発者用に提供されている欠損値などの検査関数
http://scikit-learn.org/stable/developers/utilities.html
sklern.utils には下記の関数の他，数値演算やテスト用のものがassert_all_finite
配列に無限値や欠損値が含まれていると ValueError を発生
check_array
配列の要素の型の変換，次元数や大きさの検査，無限値・欠損値の有無などをまとめて実行できる
check_X_y
check_array に加え，特徴用配列 X と教示情報用配列 y の大きさの一致を検査56

演算エラー処理

np.log(0) など有限にならない値を求めると，警告メッセージが表示される

このような演算例外に対する挙動を変更するのが np.seterr()

状況の指定：全般 all，割り算 divide，オーバーフロー over，アンダーフロー under，不正演算 invalid

動作の指定：無視 ignore，警告表示 warn/print，例外送出 raise，
関数呼び出し call，ログの記録 log

np.seterr(all=‘ignore’) とすると警告メッセージを全て抑制できる57

scipy.optimize
scipy.optimize：非線形最適化のパッケージ，機械学習ではよく利用
minimize_scolar()：有限区間上のスカラー関数の最適化
brent brent法，golden 黄金分割法
minimize()：多次元の非線形最小化
勾配を使わない最適化：勾配を求める必要はないが，代わりに目的関数の評価回数が増えるので非常に遅い
Nelder-Mead ネルダー・ミード法，Powell パウエル法
勾配を使う最適化：収束は早いが，勾配やヘシアンも与える必要
CG 共役勾配法，BFGS BFGS法，Newton-CG ニュートン法※ BFGSやニュートンは収束は早いがヘシアン計算で遅くなることも制約付きの最適化：パラメータが有限区間にある多次元最適化
L-BFGS-B，TNC 切断ニュートン法，COBYLA58

第Ⅱ部：機械学習アルゴリズムの実装

数値計算プログラミングの注意点

59

Python や NumPy の浮動小数点数
数値計算プログラミングでは，数学上の概念である無限の精度がある実数ではなく，精度が有限の浮動小数点数で計算する

浮動小数点数：IEEE 754[符号] × 2指数 × (1+仮数)
符号部 (sign)：正負の符号を表す
指数部 (exponent)：桁数を表す

仮数部 (mantissa)：有限桁の整数を表す

数値誤差：統計学では観測の過程で生じるものだが，計算機科学としては浮動小数点を用いることによる誤差もある
参考文献

伊理正夫，藤野和建『数値計算の常識』共立出版 (1985)

幸谷智紀『Python数値計算プログラミング』講談社 (2021)
60

Python と NumPy の整数型
Python の整数型は任意の桁数の数値（メモリが許す限り）を扱える
NumPy の数値は有限のbit数（int64 や uint32 など）なので，表現できる数の桁数は有限になる
大きすぎる桁数の数を代入すれば，桁あふれ (overﬂow) となる
浮動小数点に代入すれば有限精度の数になる
サンプルファイル：python_numpy_integer.ipynb







608那由多 2818阿僧祇 6403恒河沙 …

608281864034267560872252163321295376887552831379210240000000000--------------------------------------------------------------------------OverflowError
Traceback (most recent call last)<ipython-input-8-d206455b3e0e> in <module>()4 print(a)5 b = np.empty(1, dtype=np.int64)----> 6 b[:] = a
7 print(b)OverflowError: Python int too large to convert to C long61

10進数と2進数
人が収集したデータ＝10進数

計算機が処理するデータ＝2進数

整数：10進数と2進数は相互に変換可能

小数：必ずしも相互に変換可能ではない ＆ 有限桁ではさらに制限
10進数：10−n 刻みの小数
10進数

2進数：2−n 刻みの小数を表現可能

3/10

1/10

2進数

6/8

サンプルファイル：numerical_errors.ipynb

1/23

a = np.empty(1, dtype=np.float64)a[:] = 0.2
np.set_printoptions(precision=8) # default
print(a)np.set_printoptions(precision=30)
print(a)有効数字8桁では違いが分からない[ 0.2][ 0.200000000000000011102230246252]

有効数字30桁では近似だと分かる
62



 

 





まるめ誤差による失敗事例

8Dg~¦8~oUN¦5[NU¦~[NgDr¦DgR¦fD¦fU¦8Dn~¦Uv¦grr¦y~¦DNp¦D¦?NR¦fUy¦fUU¦g¦D¦DyaU¦aDU¦fhX¦~X¦UNUy¦~¦v~U¦!UNDU¦fU¦
fh_¦h¦RhUNr¦~~h~yDr¦~¦gvU¦UD~rDgya¦fU¦.DUrg¦RDD¦fgNf¦
hyRgNDUR¦D¦ UNUy¦fg_¦DXU¦f~¦RUUvgyUR¦fD¦fU¦DyaU¦aDU¦
~rR¦fg_¦ UNUy¦D_U¦DL~¦f~¦~X¦N~ygy~¦U¦?UNg\NDrr¦DXU¦DL~¦f~¦fU¦hyDNNDU¦hvU¦NDrNrDg~y¦LUN~vU¦[NhUyr¦rDaU¦~¦NDU¦fU¦DRD¦~¦r~~p¦gy¦fU¦~ya¦rDNU¦]¦fU¦DaU ¦#~yU Uyr¦fU¦Uv¦Zhr¦~¦DNp¦DyR¦hyUNU¦fU¦?NR ¦?UU¦)aU¦¦]¦D¦UUUyDg~y¦~X¦Dy¦gyN~UN¦DyaU¦aDU¦NDrNrDg~y ¦UyRg¦..¦f~¦fU¦yvLU¦~X¦vUU¦fU¦DyaU¦aDU¦fh_¦^~v¦fU¦NUyU¦~X¦fU¦
DaU¦D_U¦N~ygy~¦~UDg~y ¦Patriot Missile Defense:

Software Problem Led to System Failure at Dhahran, Saudi Arabia

パトリオットシステムが，1991
年に，サウジアラビアで，イランのミサイルの迎撃に失敗した(ibV¦¦/zOVOs¦$EsOsEVS¦<EzbV¦*EV¦

時間のカウンタを24bit精度の浮動小数点に変換

 BFOq¦Qj{¦£¦6{t ¦=HcV¦+FVT
9k{¦Y¦"VFw¦9PVT

1kmtV¦6jTV¦=HdV¦+FV¦ AVJ¦{¡
4¦>IV¦,KV¢&}kV¦"VEw:¤¥

>E|V¦
,EV¦
VF¦

2jjtV¦

時間が経つと誤差が蓄積して，時間の精度を維持できない
https://www.gao.gov/assets/220/215614.pdf-7 03C'%¦ ;G¦3lluW¦@`GW¦;MuWx¦;GeW¦
63数値誤差
量 x の測定値・浮動小数点による表現値 a に誤差 Δa＞0 がある
x ∈ (a - Δa, a + Δa)

量 y についても同様に観測値・表現値 b に誤差 Δb＞0 がある
x と y についての演算で z を計算，z の表現値は c その誤差 Δc

加減算：絶対誤差が和になる

z=x±y<latexit sha1_base64="SGOWxkdmQBzWW2mcuVOwMcAKJL0=">AAADInichVI7S8RAEP6M7/epjSBI8DgRhGNPBB8giFpop3eeHqgcSVw1mEtCsneoh52VjaWFjQoWaqn/wEawErEQ/ANiqWCj4CSXwxfqhM1+O7vfzLezo9qG7grG7kqk0rLyisqq6prauvqGxlBT84xrZR2NJzXLsJyUqrjc0E2eFLoweMp2uJJRDT6rro56+7M57ri6ZU6LdZsvZJRlU1/SNUWQKx1qnx/jhlBkTR6SA6jI3UWopkNhFmW+yT9BLABhBDZphV4wj0VY0JBFBhwmBGEDClz65hADg02+BeTJ5xDS/X2OTdQQN0unOJ1QyLtK/2VazQVek9ZeTNdna5TFoOEQU0aE3bJj9sQu2Sl7YK+/xsr7MTwt6zSrBS63043brYmXf1kZmgVWPlh/ahZYQr+vVSfttu/xbqEV+LmN3afEYDyS72SH7JH0H7A7dkE3MHPP2tEUj+9R9Ajknbf4frzij0wmVcGLvBbU1fWrulbIQqOoWiY9lv8ig4QTmMYEUp+8v9++GKFYNa/yrvdi1CCx7+3wEyR7ogPR2FRveHgk6JQqtKEDXdQOfRjGOCaRpIRbOMEZzqUj6Uq6lm4KR6WSgNOCLybdvwPTL7Wf</latexit><latexit sha1_base64="1Djs3zVSVgu2RBcfsWuCQvFLAZk=">AAADDXicSyrIySwuMTC4ycjEzMLKxs7BycXNw8vHLyAoFFacX1qUnBqanJ+TXxSRlFicmpOZlxpaklmSkxpRUJSamJuUkxqelO0Mkg8vSy0qzszPCympLEiNzU1Mz8tMy0xOLAEKRVcp2CpUKMQU5CpUxgsoG+gZgIECJsMQylBmgIKAfIHvDDEMKQz5DMkMpQy5DKkMeQwlQHYOQyJDMRBGMxgyGDAUAMViGaqBYkVAViZYPpWhloELqLcUqCoVqCIRKJoNJNOBvGioaB6QDzKzGKw7GWhLDhAXAXUqMKgaXDVYafDZ4ITBaoOXBn9wmlUNNgPklkognQTRm1oQz98lEfydoK5cIF3CkIHQhdfNJQxpDBZgt2YC3V4AFgH5Ihmiv6xq+udgqyDVajWDRQavge5faHDT4DDQB3llX5KXBqYGzQaarsqg0P83aEEQGx6b8oChADK5AhquxeBQrYDYAsQwVysA3ZMPjhErIDuYIYTBkyECSRS372EmwEINFPLFoBgDJhBD9OSAyQg10rPUMww0UXZwgqYUDgZpBiUGDWByMGdwYPBgCGAIBUfoVIZ5DPOZJjBtY9rDtA+ilIkRqkeYAQUwHQcAhnqu0Q==</latexit>

乗除算：相対誤差が和になる

z = xy or z = x_y<latexit sha1_base64="+MkFWgzxlm/ZvlF4U+3WIuYdh08=">AAADHnichVLPSxtBFP7c+qtqNdZLaS+DIdJTnIigDRSkvdSbJkYDKmF3nejiZnfZnYT8QPDcS4899FKFCtKD/g8e9CKeBP0TpEcLXhR8O9mgrahvmZ0338z33jdvnuHZViA5P2/TXrR3dHZ1v+zp7XvVPxAbfD0fuGXfFDnTtV0/b+iBsC1H5KQlbZH3fKGXDFssGOufw/2FivADy3XmZM0TyyV91bGKlqlLggqxt3X2kVVrbEmKqmww12cbTEGjtUIszpNcGXvopCInjshm3NgVlrACFybKKEHAgSTfho6AvkWkwOERtowGYT55ltoX2EAPcct0StAJndB1+q/SajFCHVqHMQPFNimLTcMnJkOCn/JdfsmP+G9+wa8fjdVQMUItNZqNJld4hYGvb7JXz7JKNEus3bGe1CxRxKTSapF2TyHhLcwmv1L/fplNZxKNEb7N/5D+LX7OD+gGTuWv+WtWZH5Q9ATYt5vMz0znE5kcqkIYuRrVNVBVrTaz0GipZqTHVS+SJj+LOUwjfw99/PatCK2qhZUPwhejBkn93w4PndxY8kMyNTsen/oUdUo33mEY76kdJjCFL5hBjhJuYgd72Ne2tEPtWDtpHtXaIs4Q/jHt7BYURbSw</latexit>

c = a+ bÛ c Û Û aÛ Û bÛ
Û Û=Û Û+Û ÛÛ c Û Û a Û Û b ÛÛ Û Û Û Û Û<latexit sha1_base64="GcjZ2q7fjVLo5WqcOfK1M6ItlBI=">AAADZ3ichVLLahRBFD2Z1hjjI62CBNy0DiOCMNRIIBoQgrrQXTKTMQPpMFSXNZMiNd1Nd82Q2OkfcOPShRsVXIif4UY/wEU+QV24iBAFBW/39PjKw2qq69atc+49det6oVaxYWx7rGQdOTp+bOL45ImTp05P2WfO3o+DfiRkUwQ6iFoej6VWvmwaZbRshZHkPU/LZW/9dna+PJBRrAJ/yWyGcrXHu77qKMENudq2cLXsmC23E3GRuHekNtwRaSJSN1LdNbPl3HT2QfA04b8QV/dDeGnijRBtu8yqLB/OXqNWGGUUYyGwd+HiAQII9NGDhA9DtgZHTN8KamAIybeKhHwRWSo/l0gxSdw+oSQhOHnX6d+l3Urh9WmfxYxztqAsmmZETAcV9p69YjvsLXvNPrDvB8ZK8hiZlk1avSFXhu2pR9ON3f+yerQarP1mHarZoIPruVZF2sPck91CDPmDh092GnP1SnKZvWCfSP9zts3e0A38wRfxclHWn1L0CpzHP+rP6uOHZPKpClnkjaKucV7VjWEWmiPVDukJ8heZI7uBJdxD6w/vwbcfRRhVLat8nL0YNUjt33bYazSvVW9Ua4sz5flbRadM4AIu4Qq1wyzmcRcLaFLCd/iMr/hW+mjZ1nlreggtjRWcc/hrWBd/Am0X0Z4=</latexit>

64

桁落ち
大きさがほぼ同じ数

誤差部分
差を計算

有効桁数は小さくなる

有効桁数

桁落ち
有効桁数に対して誤差部分の影響が相対的に拡大する
大きさがほぼ同じの大きな数の差は誤差を大きくする
65

誤差に特に注意が必要な場合
加減算：絶対誤差が和になる

大きな数の加減算で，解の絶対値が小さい（同じくらいの数の差）解に対して，絶対誤差の影響が相対的に大きい
桁落ちを生じやすい
大きな数と小さな数の加減算
小さい方の値に対して，解の誤差の影響が相対的に大きい

乗除算：相対誤差が和になる

小さな数と大きな数の乗除算
相対誤差の大きくなりやすい小さな数が，解の相対誤差を大き
くする

66

桁落ち
大きな数同士の差による桁落ちの例
本来は負にならない（二乗平均）­（平均の二乗）が負になる
サンプルファイル：numerical_errors.ipynb
n = 1000
x = 1000 + np.random.randn(n) * 10**-5
print(np.mean(x**2) - np.mean(x)**2)

1000 + 小さな乱数-1.16415321827e-10

大きな絶対値を引く
誤差に対する有効桁数を相対的に大きくして対処
y = x - 1000
print(np.mean(y**2) - np.mean(y)**2)

大きな絶対値を引く

9.50575571567e-11
67大きな数と小さな数の和大きな数と小さな数の加減算

大きな数の絶対誤差＞小さな数

小さな数の総和：最後の方は大きな数と小さな数の和は，大きな数にとって誤差となる
サンプルファイル：numerical_errors.ipynb
k = 5
a = 0
for i in range(10**k):a += 10**-k
print(1 - a)

1/10000 を 10000回 足す

1.91624494050302e-12
a = 0.
for i1 in range(10):b = 0.
for i2 in range(10):skip…
a += bprint(1 - a)-6.661338147750939e-16

小さな数同士を先に足して大きな数にする

分けて計算すれば影響は小さい
68

桁あふれ
シグモイド関数など指数を使う関数は容易に桁あふれを起こす
1
sig( ) =1 + exp( )<latexit sha1_base64="mEq1Y5XdXXsWMKhkVuJDknNUEjg=">AAADJ3ichVLPTxNBFP5YEBFRKlxMvExsSkqMzSyaSAgmJF7gBi2FJpQ0u+u0Ttjubna2zcKmVw9cOHrwIiYcCF71H/DijRNGjh6NR0z0oIlvt9v4g1DeZnbefDPfe9+8eaZnSxVwfjqgDQ5dGb46cm30+tiNm+OZWxNrym35lihbru36FdNQwpaOKAcysEXF84XRNG2xbm49iffX28JX0nVWg21PbDaNhiPr0jICgmoZVlWykQ+n2WNWrfuGFemdSGf3WFWEXv5+ON2pZbK8wBNj5x09dbJIbdnN/EAVT+HCQgtNCDgIyLdhQNG3AR0cHmGbiAjzyZPJvkAHo8Rt0SlBJwxCt+jfoNVGijq0jmOqhG1RFpuGT0yGHD/hh/yMf+BH/Av/eWGsKIkRa9mm2exyhVcb371d+n4pq0lzgGd/WH01B6hjNtEqSbuXIPEtrC6/vfPirDRXzEVT/DX/Svr3+Sl/Tzdw2t+sgxVRfEnRc2B7v4qvisN9MjlUhThymNZVJVUNMdPNQ6Onm5EiN3mTOfJLWMUSKn+hF9+/F6FXt7j2Kn4zahH9/4Y476zNFPQHBX3lYXZhPm2WEdzBXeSpIx5hAYtYRpkyPscbvMU77VA71j5qn7pHtYGUM4l/TPv8G9JCtxI=</latexit>

サンプルファイル：numerical_errors.ipynb
np.set_printoptions(precision=30)
for x in range(20, 50, 5):
print(x, sig(x))
print(sig(50) == 1.0)
20 0.9999999979388463
25 0.999999999986112
30 0.9999999999999065
35 0.9999999999999993
40 1.045 1.0
本来は 1 にならないはずだが，1 になる
True

x の入力の範囲をクリップしてしまう
x = np.clip(x, -34.538776394910684, 34.538776394910684)
69小さな数の割り算
大きな数 / 小さな数

=

特異値分解
1 =
<latexit sha1_base64="nFoMInFqVsaXopARRVXd2BrPnxU=">AAADW3ichVLPaxNBFP6a1bbG2qaKIHhZDCk9pGE2CpWCUPWitza/GsjGsLuZpEM3u8vOJLSG/gNePHrwooIH8T/w6sV/wEMPgtfqsYIeFHy72RBtafuW2fnmm/nee/Pm2YErpGLsYCqlXbg4PTN7KX157sr8Qmbxak36/dDhVcd3/bBuW5K7wuNVJZTL60HIrZ7t8i1752G0vzXgoRS+V1F7AW/2rK4nOsKxFFGtTMW0O/f1ezpNVTNvtoXVbZgu6dtWy8jrY1gk2PaVdCbUMFgx9ifLoGnmIy+1J2allcmyAotNPwmMBGSR2Iaf+QUTbfhw0EcPHB4UYRcWJH0NGGAIiGtiSFxISMT7HPtIk7ZPpzidsIjdoX+XVo2E9Wgd+ZSx2qEoLo2QlDpy7DN7x47YJ/aeHbLfp/oaxj6iXPZotkdaHrQWnt0o/zxX1aNZYXuiOjNnhQ7uxrkKyj2ImegWzkg/ePriqLxWyg2X2Bv2nfJ/zQ7YR7qBN/jhvN3kpZfkPQf9+Z/Sq9L0GZE8qkLkeTepq4yruoviKA6Ncd46ZeTHb7JGuIwKHqP+D3v6/ccexnWLai+jN6MWMY43xElQKxaM2wVj8052/UHSLLO4iVtYpo5YxToeYQNVivgBX3GIb6kvmqaltbnR0dRUormG/0y7/hci1MdZ</latexit>

逆行列<latexit sha1_base64="cCFqwcrrFVX3iNRnLCxb0YDRn7A=">AAADaHichVJNb9NAEH2NoZTyUdMeAHGxGgVxCGEdkFpVQipwgVubr0bKppHtbMKqjm3Zm6gl6h/gwpEDF0DigPgZHOAPcOhPAMSpiHIAibHjqIGqZaz1zr7d92Z2duzAlZFibG8qo506PX1m5uzsufMXLs7pl+Zrkd8PHVF1fNcP67YVCVd6oqqkckU9CIXVs12xYW89iPc3BiKMpO9V1E4gmj2r68mOdCxFUEu3uN25t8mlNzDuGuTXeJ63pdVtmLe4SzJtq2XmjcNFMW/wtq8iZxIcBjfN3UkgaPJ8rFbd5JWWnmUFlphx1DFTJ4vU1nz9Jzja8OGgjx4EPCjyXViI6GvABENAWBNDwkLyZLIvsItZ4vbplKATFqFb9O/SqpGiHq1jzShhOxTFpRES00COfWJv2T77yN6xz+zXsVrDRCPOZYdme8QVQWvu6ZXywX9ZPZoVHh+yTsxZoYPlJFdJuQcJEt/CGfEHT57vl1dKueF19pp9pfxfsT32nm7gDb47b9ZF6QWp52A8+116WZo+IZJHVYiVt9O6RklVt1EcxaExztugjPzkTVbIL6OCR6hPoMfff6wwrltc+yh+M2oR89+GOOrUigXzdsFcv5NdvZ82ywyuYRE3qCOWsIqHWEOVIn7AN/zAQeaLpmuXtaujo5mplLOAv0xb/APUksti</latexit>

小さな数の相対誤差が結果の相対誤差

diag[ 1 , 2 , … ,diag[1

1, 1

2, … , 1

1,

]

1, 1

]

割り算

逆行列の計算では固有値の除算により相対誤差が大きくなりやすい↓
小さな固有値は無視して行列のランクを小さくして回避
条件数 (conditional number)

1<latexit sha1_base64="O8fC05K6vrEwp8/NNN5TItV7jQk=">AAADHnichVJNS9xAGH6Mn7W1rnqo4GXpssVDWSYqKJ7UXvSmu64uGF2SdFYHs0lIsou6+gd66bEHL37gQWx/hSBepKeC/gTxqFAPCn0zyWJbUSdk5pln5nnnmXdew7WEHzB20aA0NjW3tLa9an/9puNtZ6Kre853Kp7J86ZjOV7B0H1uCZvnAxFYvOB6XC8bFp83Vj+F6/NV7vnCsWeDdZcvlvVlW5SEqQdEFRPvtE3NKI1T/zFCS5qwq9pmMZFiGSZb8jFQY5BC3KadxC00fIYDExWUwWEjIGxBh0/fAlQwuMQtokacR0jIdY4ttJO2Qrs47dCJXaV+mWYLMWvTPIzpS7VJp1j0e6RMIs1+sUN2zU7ZEbtkd0/GqskYoZd1Go1Iy91i55fe3O8XVWUaA6w8qJ71HKCEEelVkHdXMuEtzEhf3fh2nRvNpmsf2B67Iv+77IId0w3s6o15MMOz2xQ9jeTX++xOtuWZk2zKQhh5Lc6rL7O6hoHoHPrrvpPkyJFvMko4h1lMofAX+/T96xHqeQtz74dvRiWi/l8Qj8HcQEYdzKgzQ6mxibhY2tCH9+inihjGGCYxjbzMzD6+44eyo5woZ8rPaKvSEGt68E9Tzv8ARie1Pg==</latexit>

連立方程式 Ax = y で x の誤差は A や y の誤差のたかだか条件数倍
np.linalg.cond ：いろいろなノルムの条件数

70数値誤差への対策
誤差があることを前提とした実装
数値誤差のために，非負などの条件が成立しない場合がある

例：分散の公式（二乗平均 ­ 平均の二乗）が負になりうる
非負を保証するように max 関数などを用いる



浮動小数点の計算では等号演算子は使わない

代入だけでも，値が変わることがあり，成立しない

np.isclose(a, b) という微少な変動を許した等号判定関数がある71数値誤差への対策
明示的に比較演算をしない場合でも，暗黙的な比較がある場合加算を繰り返すことの誤差によって，大小比較の結果が変わるた
め，指定した終了値で終わらない場合がある




この場合は指定した区間を区切る  で対応できる




72

桁あふれ対策
最後の計算結果の最大化するだけなら，対数などをとる

桁あふれが生じそうな入力に対しては，定数などで置き換える
Π記号の計算で対数をとってから最後に指数関数を適用する
sp.misc.logsumexp() などの専用関数もある

途中の結果が大きく，最後に小さくなる場合
例：ユークリッド距離

大きな  で正規化しておく


73

第Ⅲ部：関連パッケージと情報源

関連パッケージ
Jupyter Notebook，scikit-statmodels，最適化など

74

Jupyter Notebook/Lab
便利な対話的な実行環境

コマンドの補完（tabキー），ヒストリの管理など

notebook は，計算結果をそのままメモとして残すことができる
ipypararrel で複数のマシン上でのジョブ管理ができる

ドキュメントの参照

オブジェクト名 + “?” で，そのオブジェクトの説明を参照できる

マジックコマンド：% (単一行)や %% (複数行) で始まるコマンド群%quickref：コマンドのリファレンス%cd，%ls：shell としてのコマンド%timeit：関数の実行時間の計測%run：ファイルの実行%pycat：カラーリング付きのソース表示75

Jupyter Notebook
分析コードとメモを一つにまとめておける

Jupyter は ipython と独立して，Python の他，R / Julia など
様々なスクリプト言語が利用できるようになった．

起動：コマンドラインでは “jupyter-notebook”

ヘルプの keyboard shortcuts を最初に読むとよい

メモの記述：Markdown 記法と，LaTeX の数式記法が利用可能
カスタマイズ：“ipython profile create” を実行すると，“ipython

profile locate” で表示される場所に“ipython_default” が作成され，カスタマイズ可能にその中の startup ディレクトリに起動ファイル .ipy を設置できる

分析結果の保存：notebook は .ipynb ファイルに保存できる
ノートブックの公開：GitHub や nbviewer https://nbviewer.ipython.org で公開できる76

IPython console for SymPy 1.0 (Python 3.5.2-6
ython)

SymPy

These commands were executed:>>> from __future__ import division>>> from sympy import *>>> x, y, z, t = symbols('x y z t')>>> k, m, n = symbols('k m n', integer=True)>>> f, g, h = symbols('f g h', cls=Function)
サンプルファイル：sympy_demo.ipynb>>> init_printing()

sympy.init_session() 実行後に行う

数式処理を行うためのパッケージ

最初に sympy.init_session() を実行
しておくと便利x, y, z, t を実数変数，k, m, n を整数変数と認識するようになる

主な機能

微分：diﬀ，積分：integrate，
極限：limit

展開：expand，単純化：
simplify，テイラー展開：series
方程式の解：solve

Documentation can be found at http://docs.sym

微分In [3]: diff((x ** 2 + log(x)) / x, x)Out[3]:

1
1
1 2
2x
+
−
(x + log (x))
x(x ) x2

積分
In [4]: integrate(x ** 3 + sin(x) ** 2, x)Out[4]:x4
x1+ − sin (x) cos (x)
4
2 2

展開
In [5]: expand((x + 1)**2)
Out[5]:

x2 + 2x + 1

77

scikit-statmodels
scikit-learn と同様に使える，scikit-learn に対する特色はARMA などの時系列分析モデル

R言語のように線形モデルを記述する記法を備える

Notebook

多重検定など検定周辺の機能

2014/11/03 6:33

サンプルファイル：statsmodels_demo.ipynb
http://statsmodels.sourceforge.net/stable/index.html のサンプルパッケージの読み込み
In [1]: import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
データの読み込みIn [2]: url = 'http://vincentarelbundock.github.io/Rdatasets/csv/HistData/Guerry.csv'
dat = pd.read_csv(url)回帰モデルのあてはめ
In [3]: results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()
結果の表示
In [4]: print results.summary()
OLS Regression Results==============================================================================
Dep. Variable:
LotteryR-squared:
0.348
Model:
OLS
Adj. R-squared:
0.333

78

Tensorflow による自動微分

自動微分

In [1]: from __future__ import (print_function,
division,
absolute_import,
unicode_literals)

サンプルファイル：tensorflow_auto_diﬀ.ipynb
Tensorflow 関連のパッケージの読み込み
In [2]: import tensorflow as tf
import numpy as np

TensorFlowなどの深層学習やAutograd などのパッケージが提供する

数値微分とは異なり，解析的に微分し，数値を代入した結果を得る
他にも数式を直接的に定義可能
非線形最適化など各種の数理計画法で結果を利用可能

スカラー変数のシンボルの定義
In [3]: x = tf.placeholder(tf.float64) # double
関数 f (x) = sin(x) を定義
In [4]: y = tf.sin(x)変数 x に，値 π を代入して関数 f (x) 評価する：sin(π) = 0
In [5]: feed_dict = {x: np.pi}
with tf.Session() as sess:
sess.run(tf.global_variables_initializer())
print(sess.run([y], feed_dict))[1.2246467991473532e-16]
In [6]: dy = tf.gradients(y, x)関数 f (x) を微分して，値 π で評価する：sin′ (π) = cos(π) = −1.0
In [7]: feed_dict = {x: np.pi}
with tf.Session() as sess:
sess.run(tf.global_variables_initializer())
print(sess.run(dy, feed_dict))[-1.0]

79

ベイズ推定Python でのベイズ推定では PyMC と PyStan が著名

どちらも，生成モデルを定義すると MCMC により，データを与えたときの変数の事後分布を推定する

PyMC3 https://docs.pymc.io/
Python でのベイズ推定パッケージとしては開発が長く続いている
Python のオブジェクトを使ってモデルを構築できる

Tensorﬂow Probability を基盤にしたPyMC4は開発頓挫
PyStan https://pystan.readthedocs.org
専用の記述言語でモデルを構築し，それを Python からでも R からでも利用できる80

超パラメータ調整
機械学習手法は適切な超パラメータの設定の必要

scikit-learn には GridSearchCV や RandomizedSearchCV がある
AutoML https://www.ml4aad.org/automl/
Optuna https://www.preferred.jp/en/projects/optuna/
性能指標の分布予測を利用した高効率な超パラメータの探索※ より多くの手法や超パラメータの中から選択するには，それに応じた十分なデータ数が必要なことには注意
サンプルファイル：auto_sklearn_demo.py
http://www.ml4aad.org/automl/
auto-sklearn/












81

最適化関連パッケージ
scipy その他の最適化関連
scipy.optimize.root 関数値を0にする入力値である根を求める
scipy.optimize.linprog 線形計画法
cvxopt http://cvxopt.org
機械学習では必要になる凸最適化に特化

凸最適化で著名な Stephen Boyd のチームが作成

TensorFlow https://www.tensorflow.org/
自動微分と確率的勾配降下による最適化

計算モデルの可視化，分散環境，ユーザコミュニティが大規模

autograd https://github.com/HIPS/autograd

TensorFlow のようにシンボルを定義せず，Python 関数そのままで自動微分ができる

82機械学習の総合パッケージ
scikit-learn https://scikit-learn.org/
よく利用され，安定している
contrib https://github.com/scikit-learn-contrib などの拡張もMLlib https://spark.apache.org/mllib/
大規模並列データ処理環境Spark上で動作する

分類・相関・クラスタリングに加え，頻出パターンマイニングもOrange https://orange.biolab.si
分析モデルを GUI を用いて対話的に構築できる統合環境
頻出パターンマイニングなどの拡張プラグインも

Vowpal Wabbit https://github.com/VowpalWabbit/vowpal_wabbit
オンライン最適化を用いた大規模・高速対応

83

深層学習

Pythonから利用できる深層学習フレームワークは多数
PyTorch http://pytorch.org/
動的に生成する計算グラフに強い
Keras https://keras.io/
TensorFlow を演算基盤に，モデルを簡易に定義可能※ 他に MXNet https://mxnet.incubator.apache.org/ や Caﬀe
https://caﬀe.berkeleyvision.org/

84

高速化
numexpr https://github.com/pydata/numexpr
要素ごとの演算をキャッシュミスが起きないようにする
bottleneck http://berkeleyanalytics.com/bottleneck/
欠損値 np.nan がある場合に高速化
Cython http://cython.org
制限はあるが Python をネイティブコードに変換できる
型宣言をコードに追加するとさらに高速化できる
高速な実行環境
PyPy：文法を制限した Python を実行
Numba：NumPy用のコンパイラ
PySton：clang などの仮想マシン LLVM 上で実行

85

並列実行環境joblib https://pythonhosted.org/joblib/jobの並列パイプライン実行

scikit-learn で利用されており，パッケージ内に含まれている
実際のバックエンドに dask や ipyparallel を利用可能

dask https://dask.pydata.org
Numpy や Pandas の操作を並列化できる
バックエンドにipyparallel を利用可能

ipyparallel https://ipyparallel.readthedocs.io
クラスタノード間の通信にいろいろな手段が使える

ローカル, SSH, MPI, PBS (SGE, Torque, LSF), Windows HPC
SSHトンネリングで，クラウド上でjobをセキュアに実行可能
jupyter notebook から操作でき，便利な %コマンド がある
86

その他のパッケージ
matplotlib https://matplotlib.org
グラフの描画
Pandas https://pandas.pydata.org
R言語のデータフレームのような表データに対する多様な操作
csvkit https://csvkit.readthedocs.org
csvファイルを操作する便利なコマンド群
Sphinx https://www.sphinx-doc.org/
ReStructuredText 形式で記述するドキュメント作成ソフト
OpenCV https://opencv.org
画像処理で著名なライブラリ．Pythonインタフェースあり
NLTK https://www.nltk.org
自然言語処理のPythonパッケージ87

第Ⅲ部：関連パッケージと情報源

情報源
Pythonによるデータ分析についての資料や勉強会

88

国際会議・勉強会
SciPy Conferences https://conference.scipy.org/
科学技術計算での Python 利用に関する国際会議

US，欧州，各国の会議情報があり，スライド・ビデオを公開

PyCon JP https://www.pycon.jp
毎年9月に開催．アジア太平洋の PyCon APAC として開催されることもPython全般の国内会議．科学技術計算系の講演が増えている

PyData Tokyo https://pydatatokyo.connpass.com
データ分析でのPython利用の勉強会
Start Python Club https://startpython.connpass.com
Pythonでスタートする人たちの集い
89拙著のSciPy関連資料
拙著のPythonプログラミング関係を含む資料

http://www.kamishima.net/jp/kaisetsu/

機械学習のPythonとの出会い (Machine Learning Meets Python)
機械学習とPythonの基礎を知っている人向けの，NumPy /SciPy を使った数値計算プログラムのチュートリアルです

具体的なアルゴリズムを実装する過程を通じて，パッケージをどう使うかを知ることができるように工夫しました
HTML版の他PDF版とePub版も配布しています
Python による科学技術計算の概要

Pythonを使った科学技術計算の講演資料をまとめたものです

朱鷺の杜Wiki http://ibisforest.org/index.php?python
関連情報をいろいろ集めてまとめています90