Pythonで直積，順列，組み合わせ，階乗の計算

●全要素のイテレータが欲しい場合（import itertools）

直積（デカルト積）：itertools.product

例：

list(itertools.product([‘a’, ‘b’, ‘c’, ‘d’, ‘e’], repeat=2))

#= itertools.product([‘a’, ‘b’, ‘c’, ‘d’, ‘e’], [‘a’, ‘b’, ‘c’, ‘d’, ‘e’])

[(‘a’, ‘a’),
(‘a’, ‘b’),
(‘a’, ‘c’),
(‘a’, ‘d’),
(‘a’, ‘e’),
(‘b’, ‘a’),
(‘b’, ‘b’),
(‘b’, ‘c’),
(‘b’, ‘d’),
(‘b’, ‘e’),
(‘c’, ‘a’),
(‘c’, ‘b’),
(‘c’, ‘c’),
(‘c’, ‘d’),
(‘c’, ‘e’),
(‘d’, ‘a’),
(‘d’, ‘b’),
(‘d’, ‘c’),
(‘d’, ‘d’),
(‘d’, ‘e’),
(‘e’, ‘a’),
(‘e’, ‘b’),
(‘e’, ‘c’),
(‘e’, ‘d’),
(‘e’, ‘e’)]

順列：itertools.permutations

例：

list(itertools.permutations([‘a’, ‘b’, ‘c’, ‘d’, ‘e’], 2))

[(‘a’, ‘b’),
(‘a’, ‘c’),
(‘a’, ‘d’),
(‘a’, ‘e’),
(‘b’, ‘a’),
(‘b’, ‘c’),
(‘b’, ‘d’),
(‘b’, ‘e’),
(‘c’, ‘a’),
(‘c’, ‘b’),
(‘c’, ‘d’),
(‘c’, ‘e’),
(‘d’, ‘a’),
(‘d’, ‘b’),
(‘d’, ‘c’),
(‘d’, ‘e’),
(‘e’, ‘a’),
(‘e’, ‘b’),
(‘e’, ‘c’),
(‘e’, ‘d’)]

組み合わせ：itertools.combinations

例：

list(itertools.combinations([‘a’, ‘b’, ‘c’, ‘d’, ‘e’], 2))

[(‘a’, ‘b’),
(‘a’, ‘c’),
(‘a’, ‘d’),
(‘a’, ‘e’),
(‘b’, ‘c’),
(‘b’, ‘d’),
(‘b’, ‘e’),
(‘c’, ‘d’),
(‘c’, ‘e’),
(‘d’, ‘e’)]

itertoolsはイテレータの操作には便利だが，数値計算には向いていないし，そもそも，別に速くもない．ループ構造をシンプルにしたり，特に直積（デカルト積）を求めたりする時よく使うけど（productは色々な場面で使う），単に直積を求める場合は，内包表記で2重ループの方が速かったりする．メモリ効率が良くなる以上の効果って余り無い気がする．

数が欲しい時，「len(list(itertools.permutations(range(365), k)))」みたいに書けるが，計算コストが尋常じゃない．

●数が欲しい（数値計算がしたい）場合

（import numpy, scipy.special, scipy.misc）

数値計算をする時は，まずSciPyを漁る（SciPy⊃Numpy）．「scipy.special.perm()」「scipy.special.comb()」（他にscipy.misc.comb()とか）なんかがそう．numpyの方が有名過ぎて，ググってもnumpyに関する情報しか出てこない事もあるけど，まあその場合は大人しくnumpyを使う．

直積： np.outer or ブロードキャスト（data[:, np.newaxis] * data）

例：

import numpy as np

a = [1, 2]

b = [3, 4]

np.outer(a, b)

array([[3, 4],
[6, 8]])

例：

np.array(a)[:, np.newaxis] * b

array([[3, 4],
[6, 8]])

順列：scipy.special.perm

例：

import scipy.special as ss

ss.perm(10, 2)

90.0

例：

ss.perm(a+b, 2)

array([ 0., 2., 6., 12.])

（a=[1, 2], b=[3, 4], a + b→[1, 2, 3, 4]，リストの足し算はオブジェクトの連結）

組み合わせ：scipy.special.comb/scipy.misc.comb

例：

ss.comb(10, 2)

45.0

例：

ss.comb(a+b, 2)

array([ 0., 1., 3., 6.])

階乗：math.factorial/scipy.misc.factorial

例：

import math

[math.factorial(i) for i in (a+b)]

[1, 2, 6, 24]

例：

import scipy.misc as sm

sm.factorial(a+b)

array([ 1., 2., 6., 24.])

●一例：誕生日のパラドックス

参考：

誕生日のパラドックス-Wikipedia

問題：

クラスに誕生日が同じ人がいる確率はどれくらいだろう．

クラスの人数が何人くらいだったらいると期待されるだろう．

回答例：

クラスの人数と累積確率のグラフを描いてみる．

n人の中で同じ誕生日の人が少なくとも2人いる場合の確率は，下式で求められる．

EQUATION: Equation　ー（1）

import numpy as np
import scipy.special as ss
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

n = np.array([365]*100)
k = np.arange(100)+1
CPD = 1 - ss.perm(n, k)/365.0**k

plt.plot(CPD)
plt.xlabel('Number of People')
plt.ylabel('Probability of a pair')

with plt.xkcd():
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(CPD)
    ax.plot([23, 23], [0, 0.5], '--')
    ax.plot([0, 23], [0.5, 0.5], '--')
    ax.set_xlabel('Number of People')
    ax.set_ylabel('Probability of a pair')
    plt.annotate(
        '23 people in a room \n the probability are over 50%',
        xy=(22, 0.507), arrowprops=dict(arrowstyle='->'), xytext=(30, 0.2))

plt.show()

xkcd()で遊ぶ時は，with構文を使うのが良い．前に使った時は，戻し方が分からなくて，結局Jupyter再起動したので（関連：一日一Python：色々なグラフの描き方）．