Up

Pythonサンプルプログラム

岡本安晴「統計学を学ぶための数学入門［上］」培風館の内容に合わせて

Pythonサンプルプログラムを以下のように用意してみた。

「下巻」用のウェブサイトも用意しています。

１章 四則演算

1.1 足し算と引き算

>>> (3 + 5) + 7

15

>>> 3 + (5 + 7)

15

>>> (3 - 5) - 7

-9

>>> 3 - (5 - 7)

5

>>>

1.2 掛け算と割り算

>>> (3 * 5) * 7

105

>>> 3 * (5 * 7)

105

>>> (5 / 2) / 4

0.625

>>> 5 / (2 / 4)

10.0

>>>

２章 順列と組み合わせ

2.1　順列

順列を求めるスクリプト例をリスト2.1に示す

リスト2.1　順列を求めるスクリプト例

import itertools

results = itertools.permutations(['Dog', 'Cat', 'Bird'], 2)

for v in results:

print(v)

リスト2.1を実行すると以下の出力を得る。

('Dog', 'Cat')

('Dog', 'Bird')

('Cat', 'Dog')

('Cat', 'Bird')

('Bird', 'Dog')

('Bird', 'Cat')

組み合わせを求めるスクリプト例をリスト2.2に示す。

リスト2.2　組み合わせを求めるスクリプト例

import itertools

results = itertools.combinations(['Dog', 'Cat', 'Bird'], 2)

for v in results:

print(v)

リスト2.2を実行すると以下の結果を得る。

('Dog', 'Cat')

('Dog', 'Bird')

('Cat', 'Bird')

３章 集合と関数

集合演算のスクリプト例を書いてみた。

リスト３　集合演算

A = set([1, 2, 3])

B = set([3, 4, 5])

print('A =', A)

print('B =', B)

print('A + B =', A | B)

print('A and B =', A & B)

print('A - B =', A - B)

リスト３を実行すると以下の結果を得る。

A = {1, 2, 3}

B = {3, 4, 5}

A + B = {1, 2, 3, 4, 5}

A and B = {3}

A - B = {1, 2}

４章　数列と極限

Napier数を無限級数の和の式

で求めるスクリプトを書いてみた。

リスト４　Napier数の計算

s = 1.0

i = 1

f = 1

while True:

s0 = s

f *= i

s += 1 / f

if s0 == s:

break

i += 1

print(s)

import numpy

print(numpy.e)

実行すると以下の結果を得る。

2.7182818284590455

2.718281828459045

５章　確率その１

確率は、多数回試行を繰り返したときの期待相対度数を表す（大数の法則）と解釈される。

サイコロ投げにおいて、１の出る確率は1／6、偶数の目が出る確率は0.5である。サイコロ投げの回数を増やしていったときの「１の目」および「偶数の目」が出た回数をそれぞれの相対度数でプロットして大数の法則を確認するシミュレーションプログラムをリスト５のように作成した。

リスト５　大数の法則のシミュレーション

import numpy as np

from numpy.random import randint

import matplotlib.pyplot as plt

n_simu = 2000

n_one = 0

n_even = 0

trace_one = []

trace_even = []

for i in range(n_simu):

v = randint(low = 1, high = 6 + 1)

if v == 1:

n_one += 1

trace_one.append(n_one / (i + 1))

if v % 2 == 0:

n_even += 1

trace_even.append(n_even / (i + 1))

plt.plot(range(1, n_simu + 1), trace_even, label = 'even')

plt.plot(range(1, n_simu + 1), trace_one, label = 'one')

plt.hlines(1/6, 0.0, n_simu)

plt.hlines(0.5, 0.0, n_simu)

plt.yticks([0.0, 1/6, 0.5, 1.0], ['0', '1/6', '0.5', '1'])

plt.legend(fontsize = 14)

plt.title('Law of Large Numbers', fontsize = 16)

plt.show()

リスト５のスクリプトを実行すると、図５のグラフが描かれる。

図５　大数の法則

６章　３角関数

関数sinとcosのグラフを描くスクリプトをリスト6.1のように書いてみた。

リスト6.1　関数sinとcosのグラフ

import math

import matplotlib.pyplot as plt

x_values = []

v_sin = []

v_cos = []

x = -0.1

while True:

x_values.append(x)

y = math.sin(x)

v_sin.append(y)

y = math.cos(x)

v_cos.append(y)

x += 0.01

if x > 6.5:

break

plt.plot(x_values, v_sin, label = 'sin')

plt.plot(x_values, v_cos, label = 'cos')

plt.xticks([0.0, math.pi/2, math.pi, math.pi*3/2, math.pi*2],

['0', '$\pi$/2', '$\pi$', '3$\pi$/2', '2$\pi$'])

plt.hlines(0.0, -0.1, 6.5)

plt.vlines(0.0, -0.5, 1.1)

plt.vlines(math.pi, -1.1, 1.1)

plt.vlines(math.pi*2, -1.1, 1.1)

plt.legend(fontsize = 16)

plt.show()

リスト6.1を実行すると図6.1のグラフが描かれる。

図6.1　関数sinとcos

次のリスト6.2は関数tanを描くものである。

リスト6.2　関数tanのグラフ

import math

import matplotlib.pyplot as plt

x_values = []

v_tan = []

x = -math.pi/2 + 0.2

while True:

x_values.append(x)

v_tan.append(math.tan(x))

x += 0.02

if x >= math.pi/2 - 0.2:

break

plt.figure(figsize = (5,7))

plt.plot(x_values, v_tan, label = 'tan')

plt.hlines(0.0, -1.7, 1.7)

plt.vlines(-math.pi/2, -5, 4)

plt.vlines(0.0, -5, 5)

plt.vlines(math.pi/2, -5, 5)

plt.xticks([-math.pi/2, 0.0, math.pi/2],

['-$\pi$/2', '0', '$\pi$/2'],

fontsize = 14)

plt.legend(fontsize = 16)

plt.show()

リスト6.2を実行すると図6.2のグラフが描かれる。

図6.2　関数tanのグラフ

７章　指数関数と対数関数

指数関数とのグラフを描くスクリプトをリスト7.1のように書いてみた。

リスト7.1　指数関数とのグラフ

import math

import matplotlib.pyplot as plt

x_values = []

v_2 = []

v_0p5 = []

x = -3.5

while True:

x_values.append(x)

y = 2.0 ** x

v_2.append(y)

y = (0.5) ** x

v_0p5.append(y)

x += 0.05

if x > 3.5:

break

plt.plot(x_values, v_2, label = '2 ** x')

plt.plot(x_values, v_0p5, label = '0.5 ** x')

plt.hlines(0.0, -3.6, 3.6)

plt.vlines(0.0, 0.0, 8)

plt.legend(fontsize = 16)

plt.xlabel('x', fontsize = 14)

plt.show()

リスト7.1を実行すると図7.1のグラフが描かれる。

図7.1　指数関数とのグラフ

対数関数とのグラフを描くスクリプトをリスト7.2のように用意してみた。

リスト7.2　対数関数とのグラフ

import math

import matplotlib.pyplot as plt

x_values = []

v_log2 = []

v_log0p5 = []

x = 0.01

while True:

x_values.append(x)

y = math.log2(x)

v_log2.append(y)

y = math.log(x, 0.5)

v_log0p5.append(y)

x += 0.05

if x > 3.5:

break

plt.plot(x_values, v_log2, label = 'log2(x)')

plt.plot(x_values, v_log0p5, label = 'log(x, 0.5)')

plt.vlines(0.0, -6.8, 6.8)

plt.hlines(0.0, 0.0, 3.6)

plt.legend(fontsize = 16)

plt.xlabel('x', fontsize = 14)

plt.show()

リスト7.2を実行すると図7.2のグラフが描画される。

図7.2　対数関数とのグラフ

８章　微分と積分

関数とその導関数のグラフを描くスクリプトをリスト8.1のように書いてみた。

リスト8.1　関数とその導関数の描画

import numpy as np

import matplotlib.pyplot as plt

def f(x):

return x ** 4 - 2 * x ** 2 + 1

def df(x):

return 4 * x ** 3 -4 * x

x = np.arange(-1.5, 1.5, 0.1)

v_f = [f(t) for t in x]

v_df = [df(t) for t in x]

plt.plot(x, v_f, label = 'f(x)')

plt.plot(x, v_df, label = 'df(x)')

plt.hlines(0.0, -1.5, 1.5)

plt.vlines(-1.0, -8, 3.5)

plt.vlines(1.0, -8, 3.5)

plt.xlabel('x', fontsize = 14)

plt.legend(fontsize = 12)

plt.show()

リスト8.1を実行すると図8.1のグラフが描かれる。

図8.1　関数とその導関数

積分

を求めるプログラムをリスト8.2のように書いてみた。

リスト8.2　積分の計算

import scipy.integrate as sciint

import matplotlib.pyplot as plt

import math

v = sciint.quad(lambda x: (1 - x ** 2) ** 0.5, -1, 1)

print(v[0])

print(math.pi/2)

リスト8.2を実行すると以下の出力を得る。

1.5707963267948983

1.5707963267948966

９章　確率その２

次式の正規分布

のグラフを描くプログラムをリスト9.1のように書いてみた。

リスト9.1　標準正規分布の描画

import scipy.stats as ss

import numpy as np

import matplotlib.pyplot as plt

x = np.arange(-3, 3, 0.01)

y = ss.norm.pdf(x)

plt.plot(x, y)

plt.show()

実行すると図9.1のグラフが描かれる。

図9.1　標準正規分布

次式の２変量の正規分布

のグラフを描くプログラムをリスト9.2のように書いてみた。

リスト9.2　2変量正規分のグラフ

import scipy.stats as ss

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

import numpy as np

m_n = ss.multivariate_normal([0.0, 0.0], [[1.0, 0.9], [0.9, 1.0]])

X, Y = np.mgrid[-3.0:3.0:0.05, -3.0:3.0:0.05]

pos = np.dstack((X, Y))

fig = plt.figure()

ax = Axes3D(fig)

ax.plot_surface(X, Y, m_n.pdf(pos), cmap = plt.cm.coolwarm, rcount = 200, ccount = 200)

plt.show()

リスト9.2を実行すると図9.2のグラフが描かれる。

図9.2　2変量正規分布

図9.2のグラフは、プログラムを実行して表示されているウィンドウ内をクリックしてドラッグすると回転する。

10章　乱数

次式の正規分布

に従う乱数を10000個生成してヒストグラムを正規分布のグラフとともに描くプログラムをリスト10.1のように書いてみた。

リスト10.1　正規乱数

import scipy.stats

import matplotlib.pyplot as plt

import numpy as np

samples = scipy.stats.norm.rvs(size = 10000)

plt.hist(samples, bins = 20, density = True)

x = np.arange(-3.0, 3.0, 0.1)

y = scipy.stats.norm.pdf(x)

plt.plot(x, y, linewidth = 3)

plt.show()

リスト10.1を実行すると図10.1のグラフが描かれる。

図10.1　正規乱数

次式の２変量の正規分布

に従う乱数の生成と正規分布の等高線グラフを描くプログラムをリスト10.2のように書いてみた。

リスト10.2　2変量正規分の乱数

import scipy.stats as ss

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

import numpy as np

m_n = ss.multivariate_normal([0.0, 0.0], [[1.0, 0.9], [0.9, 1.0]])

samples = m_n.rvs(size = 500)

X, Y = np.mgrid[-4.0:4.0:0.05, -4.0:4.0:0.05]

pos = np.dstack((X, Y))

fig, ax = plt.subplots()

ax.contour(X, Y, m_n.pdf(pos))

ax.plot(samples.T[0], samples.T[1], 'o', alpha = 0.3)

plt.title(r'$\rho$ = 0.9', fontsize = 20)

plt.show()

リスト10.2をじっくすると図10.2のグラフが描かれる。

図10.2　2変量正規分布

なお、乱数（MCMCサンプリング）を用いてベイズ分析における事後分布を推定するプログラム例は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第３部「ベイズ分析」で説明している。

11章　行列と基本演算

逆行列の計算例のプログラム例をリスト11.1のように書いてみた。

リスト11.1　逆行列の計算

import numpy as np

A = np.array([[1, 1], [0, 1]])

print('A =\n', A)

invA = np.linalg.inv(A)

print('invA =\n', invA)

print('invA * A =\n', invA @ A)

print('A * invA =\n', A @ invA)

リスト11.1を実行すると以下の出力を得る。

A =

[[1 1]

[0 1]]

invA =

[[ 1. -1.]

[ 0. 1.]]

invA * A =

[[1. 0.]

[0. 1.]]

A * invA =

[[1. 0.]

[0. 1.]]

なお、行列演算のプログラム例は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第4章「行列演算とPythonスクリプト」および第２部「多変量解析」で説明している。

Up