簡単な数学計算Pythonサンプルスクリプト
岡本安晴「数学を学ぶための数学入門［下］」培風館で扱われている内容に関連する
Pythonサンプルスクリプトを用意してみた。
「上巻」用のウェブサイトも用意している。
なお、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞において、行列演算とPythonスクリプト（第４章）、確率計算とPythonスクリプト（第9章）の説明を行っている。
第１章 行列の基礎概念
行列はNumpyのndarray型を使うと便利である。
行列
は以下のようにndarray型（n-dimensional array）のオブジェクトとして生成できる。
>>> import numpy
>>> A = numpy.array([[3, 1,
2],
[5, 1, 7]])
>>> A
array([[3, 1, 2],
[5, 1,
7]])
>>>
行列の加減算およびスカラーとの積は以下のようにできる。
>>> import numpy
>>> A = numpy.array([[3, 1,
2], [5, 2, 7]])
>>> B = numpy.array([[5, 7,
1], [6, 2, 3]])
>>> A + B
array([[ 8, 8,
3],
[11, 4, 10]])
>>> A - B
array([[-2, -6, 1],
[-1, 0,
4]])
>>> C = numpy.array([[1, 2],
[3, 4]])
>>> 3 * C
array([[ 3, 6],
[ 9, 12]])
>>>
Pythonのリストとnumpyのndarray型の演算との違いは、例えば以下のようである。
>>> [1, 2] * 3 + [4, 5]
[1, 2, 1, 2, 1, 2, 4, 5]
>>> import numpy
>>> numpy.array([1, 2]) * 3 +
numpy.array([4, 5])
array([ 7, 11])
>>>
ndarray型の要素の条件の判定結果を要素とするndarray型は以下のように生成できる。
>>> import numpy as np
>>> a =
np.array([9,1,8,3,7,2,5,6,4])
>>> check_a = a % 2 == 0
>>> print(check_a)
[False False True False False True False True True]
条件の判定結果がTrueである要素は以下のように選ぶことが出来る。
>>> b = a[check_a]
>>> print(b)
[8 2 6 4]
リストでは、ndarray型のように要素ごとの条件の判定結果を同じような操作で得ることはできない。行うと、以下のようにエラーとなる。
>>> L = [9,1,8,3,7,2,5,6,4]
>>> check_L = L % 2 == 0
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: unsupported operand type(s)
for %: 'list' and 'int'
>>>
リストとndarray型についてはこのウェブサイトでも説明している。
転置行列は以下のように得ることができる。
>>> import numpy
>>> A = numpy.array([[1, 2],
[3, 4]])
>>> A.transpose()
array([[1, 3],
[2, 4]])
>>> numpy.transpose(A)
array([[1, 3],
[2, 4]])
>>> A.T
array([[1, 3],
[2, 4]])
>>>
内積は以下のように計算できる。
>>> import numpy
>>> a = numpy.array([1, 2,
3])
>>> b = numpy.array([10, 100,
1000])
>>> a * b
array([ 10,
200, 3000])
>>> a @ b.T
3210
>>> numpy.inner(a, b)
3210
>>>
単純に積の記号 * を用いると、要素ごとの積が計算されるが、行列の積記号 @ 、あるいは関数innerを用いると内積が計算される。
行列の積の演算子 @ は以下のように計算される。
>>> a.T @ b
3210
>>> a = numpy.array([a])
>>> b = numpy.array([b])
>>> a.T @ b
array([[ 10,
100, 1000],
[ 20,
200, 2000],
[ 30,
300, 3000]])
>>> a
array([[1, 2, 3]])
>>> b
array([[ 10,
100, 1000]])
>>> a @ b.T
array([[3210]])
>>>
ベクトル（１次元ndarray型）同士の@による積は内積であるが、２次元のndarray型同志の積の場合は、積の結果は行列である。array([[3210]])は１行１列の２次元ndarray型である。
行列の要素を１列に並べる演算子としてvec演算子と呼ばれているものがあるが、numpyにも行列の要素を１列に並べる演算がある。行単位に並べるものと列単位に並べるものがあるので注意が必要である。以下に例を示す。
>>> a = numpy.array([[1, 2],
[3, 4]])
>>> a
array([[1, 2],
[3, 4]])
>>> a.reshape(-1)
array([1, 2, 3, 4])
>>> a.flatten()
array([1, 2, 3, 4])
>>> a.flatten('F')
array([1, 3, 2, 4])
>>> a.flatten('C')
array([1, 2, 3, 4])
>>>
関数reshapeの引数は要素の個数を設定するが、個数が分からないときは−1を設定すると要素数に合わせた処理が行われる。
関数flattenの引数に’F’を設定すると、列単位で並べられるが、文字’C’を設定すると行単位になる。
クロネッカー積（Kronecker product）・直積（direct product）は、関数kronにより以下のように計算できる。
>>> import numpy
>>> A = numpy.array([[2, 10],
[100, 1000]])
>>> B = numpy.array([[1, 2],
[3, 4]])
>>> numpy.kron(A, B)
array([[ 2, 4, 10, 20],
[ 6, 8, 30, 40],
[
100, 200, 1000, 2000],
[
300, 400, 3000, 4000]])
>>>
トレースは以下のように求めることができる。
>>> import numpy
>>> A = numpy.array([[1, 2],
[3, 4]])
>>> A
array([[1, 2],
[3, 4]])
>>> numpy.trace(A)
5
>>>
ランク（階数）は以下のように求めることができる。
>>> from numpy.linalg import
matrix_rank
>>> import numpy
>>> A = numpy.array([[1, 0],
[0, 1], [1, 1], [1, 2]])
>>> A
array([[1, 0],
[0, 1],
[1, 1],
[1, 2]])
>>> matrix_rank(A)
2
>>>
ノルムは以下のように求めることができる。
>>> import numpy
>>> from numpy.linalg import
norm
>>> v = numpy.array([1, 2])
>>> v
array([1, 2])
>>> norm(v)
2.23606797749979
>>> A = numpy.array([[1, 2,
3], [4, 5, 6]])
>>> A
array([[1, 2, 3],
[4, 5,
6]])
>>> norm(A)
9.539392014169456
>>>
なお、行列演算については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第４章でも説明している。
第２章 逆行列・一般逆行列
逆行列は以下のように求めることができる。
>>> import numpy
>>> from numpy.linalg import inv
>>> A = numpy.array([[1, 2],
[3, 4]])
>>> A
array([[1, 2],
[3, 4]])
>>> inv(A)
array([[-2. , 1. ],
[ 1.5,
-0.5]])
>>> A @ inv(A)
array([[1.0000000e+00, 0.0000000e+00],
[8.8817842e-16, 1.0000000e+00]])
一般化逆行列・一般逆行列（generalized inverse, g inverse）は、ムーア・ペンローズ逆行列（Moore-Penrose inverse）として以下のように求められる。
>>> import numpy
>>> from numpy.linalg import
pinv
>>> A = numpy.array([[1, 2,
3], [4, 5, 6]])
>>> A
array([[1, 2, 3],
[4, 5,
6]])
>>> pinv(A)
array([[-0.94444444, 0.44444444],
[-0.11111111, 0.11111111],
[
0.72222222, -0.22222222]])
>>> A @ pinv(A) @ A
array([[1., 2., 3.],
[4., 5.,
6.]])
>>> pinv(A) @ A @ pinv(A)
array([[-0.94444444, 0.44444444],
[-0.11111111, 0.11111111],
[
0.72222222, -0.22222222]])
>>> A @ pinv(A)
array([[ 1.00000000e+00, 2.22044605e-16],
[-8.88178420e-16,
1.00000000e+00]])
>>> pinv(A) @ A
array([[ 0.83333333, 0.33333333, -0.16666667],
[
0.33333333, 0.33333333, 0.33333333],
[-0.16666667, 0.33333333, 0.83333333]])
>>>
第３章 行列式
行列式は、以下のように求めることができる。
>>> import numpy
>>> from numpy.linalg import
det
>>> A = numpy.array([[1, 1],
[-1, 1]])
>>> A
array([[ 1, 1],
[-1, 1]])
>>> det(A)
2.0
>>>
第４章 固有値と特異値
固有値は以下のように求めることができる。
>>> import numpy
>>> from numpy.linalg import
eig
>>> A = numpy.array([[5, 3],
[3, 5]])
>>> A
array([[5, 3],
[3, 5]])
>>> w, v = eig(A)
>>> w
array([8., 2.])
>>> v
array([[ 0.70710678, -0.70710678],
[
0.70710678, 0.70710678]])
>>> for i in range(2):
print('Eigen
value-{0} =\n {1}'.format(i, w[i]))
print('Eigen
vector-{0} =\n {1}'.format(i, v.T[i]))
Eigen value-0 =
8.0
Eigen vector-0 =
[0.70710678 0.70710678]
Eigen value-1 =
1.9999999999999991
Eigen vector-1 =
[-0.70710678 0.70710678]
>>> from numpy.linalg import
eigh
>>> wh, vh = eigh(A)
>>> wh
array([2., 8.])
>>> vh
array([[-0.70710678, 0.70710678],
[
0.70710678, 0.70710678]])
>>> for i in range(2):
print('Eigen
value-{0} =\n {1}'.format(i, wh[i]))
print('Eigen
vector-{0} =\n {1}'.format(i, vh.T[i]))
Eigen value-0 =
2.0
Eigen vector-0 =
[-0.70710678 0.70710678]
Eigen value-1 =
8.0
Eigen vector-1 =
[0.70710678 0.70710678]
>>> B = numpy.array([[1, 2],
[3, 4]])
>>> B
array([[1, 2],
[3, 4]])
>>> eig(B)
(array([-0.37228132, 5.37228132]), array([[-0.82456484,
-0.41597356],
[
0.56576746, -0.90937671]]))
>>>
特異値分解（singular value decomposition: SVD）は以下のように計算できる。
>>> import numpy
>>> from numpy.linalg import
svd
>>> A = numpy.array([[2, 1,
0], [0, 1, 2]])
>>> A
array([[2, 1, 0],
[0, 1,
2]])
>>> u, s, vh = svd(A)
>>> u
array([[-0.70710678, -0.70710678],
[-0.70710678, 0.70710678]])
>>> s
array([2.44948974, 2. ])
>>> vh
array([[-5.77350269e-01,
-5.77350269e-01, -5.77350269e-01],
[-7.07106781e-01, 2.22044605e-16, 7.07106781e-01],
[
4.08248290e-01, -8.16496581e-01,
4.08248290e-01]])
>>> u @ numpy.diag(s[0:2]) @
vh[0:2]
array([[ 2.00000000e+00, 1.00000000e+00, 2.22044605e-16],
[-5.55111512e-16,
1.00000000e+00,
2.00000000e+00]])
>>>
第５章 行列と線形写像
線形写像を表す行列の行列式は写像の符号付き拡大率を表す。
行列
に対して変換行列
による変換
において確認するプログラムを以下のように作成してみた。変換は右手系を左手系に変換するものであるので、変換行列の行列式は負の値である。
>>> import numpy
>>> from numpy.linalg import
det
>>> A = numpy.array([[2, 1], [1,
2]])
>>> det(A)
2.9999999999999996
>>> T = numpy.array([[2, 3],
[3, 2]])
>>> T @ A
array([[7, 8],
[8, 7]])
>>> det(T @ A)
-14.999999999999993
>>> det(T @ A) / det(A)
-4.999999999999998
>>> det(T)
-5.000000000000001
>>>
第６章 微分
微分係数をtensorflow 2を使って求めてみた。関数は
とした。このとき
である。
スクリプトは
import tensorflow as tf
import numpy as np
def f(x):
return x ** 3
for i in np.arange(0.0, 3.1, 1.0):
x = tf.Variable(i)
with tf.GradientTape() as
tape:
y =
f(x)
df_dx = tape.gradient(y, x)
print('df/dx = ',
df_dx.numpy(), ' 3 * X**2 =
', 3 * i**2)
とした。実行すると
df/dx = 0.0 3 * X**2 = 0.0
df/dx = 3.0 3 * X**2 = 3.0
df/dx = 12.0 3 * X**2 = 12.0
df/dx = 27.0 3 * X**2 = 27.0
と出力される。
第７章 ベクトル・行列関数の微分表記
偏微分係数をtensorflow 2を使って求めてみた。
関数は
とした。このとき
である。
計算例として以下のスクリプトを用意した。
import tensorflow as tf
import numpy as np
def func(X):
return X[0] * (X[1] **
2)
for i1 in np.arange(0.0, 2.1, 1.0):
for i2 in np.arange(0.0,
2.1, 1.0):
x =
tf.Variable([i1, i2])
with
tf.GradientTape() as tape:
y = func(x)
df_dx = tape.gradient(y, x)
print(df_dx.numpy(), '
i2**2 = ', i2**2, ' 2*i1*i2
= ', 2*i1*i2)
実行すると以下の出力を得る。
[0. 0.] i2**2 = 0.0 2*i1*i2 = 0.0
[1. 0.] i2**2 = 1.0 2*i1*i2 = 0.0
[4. 0.] i2**2 = 4.0 2*i1*i2 = 0.0
[0. 0.] i2**2 = 0.0 2*i1*i2 = 0.0
[1. 2.] i2**2 = 1.0 2*i1*i2 = 2.0
[4. 4.] i2**2 = 4.0 2*i1*i2 = 4.0
[0. 0.] i2**2 = 0.0 2*i1*i2 = 0.0
[1. 4.] i2**2 = 1.0 2*i1*i2 = 4.0
[4. 8.] i2**2 = 4.0 2*i1*i2 = 8.0
第８章 積分
次の積分
の計算を行うスクリプトを以下のように書いてみた。
from scipy.integrate import quad
import math
def f(x):
return (1.0 - x ** 2) ** 0.5
v = quad(f, -1.0, 1.0)
print(v)
print(math.pi/2)
実行すると以下の結果となる。
(1.5707963267948983,
1.0002354500215915e-09)
1.5707963267948966
>>>
次の２重積分
の積分を求めるスクリプトを以下のように書いてみた。
from scipy.integrate import dblquad
import math
def f(x, y):
return 1.0
def lower_f(x):
return -(1.0 - x ** 2) **
0.5
def upper_f(x):
return (1.0 - x ** 2) ** 0.5
v = dblquad(f, -1, 1.0, lower_f,
upper_f)
print(v)
print(math.pi)
実行すると以下の結果を得る。
(3.1415926535897967,
2.000470900043183e-09)
3.141592653589793
>>>
第９章 確率
サンプルスクリプトを＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞第９章「確率計算とPythonスクリプト」に用意している。