確率論のσ加法族（ほんとは有限だけど）の自動生成を Sympy で書いてみた！

## はじめに

[「確率空間」を学生に教えるとき](https://ktysd.github.io/lec/kalman/)，彼我ともに，最も難しく感じるのは，σ加法族（完全加法族，σ代数）のところです．
そこで，有限集合に話を限って，

+ σ加法族を自分達で作らせる（本当は有限加法族だけど）

という教え方をするのですが，１つ問題がありまして，解答例を作るのが面倒なのです．
暗算しながら板書していくと，「それだと足りなくないですか」とか言われて「あっ」みたいなことになる．教室総出の人海戦術が始まり，「・・・これでいいんだっけ，まだ足んないか」みたいになって，しばし中断．そういう場面も有意義だとは思うのですが，確証ある正解は，やはり手元に欲しいものです．
というわけで，授業中の答え合わせを念頭に，

+ σ加法族（ほんとは有限だけど）の自動生成を，Python3+Sympy で書いてみた！

というのが，この記事の内容です．

確率空間なんて興味ない人のほうが多い気がするので，最低限の基本事項をレビューしながら話を進めます．

## 基本事項 ー 確率空間と`FiniteSet`

確率空間というのは，$(\Omega,{\mathcal F},P)$という三つ組のことです．$\Omega$を見本空間（または標本空間），${\mathcal F}$をσ加法族（または完全加法族，σ代数）といいます．$P$は確率です（この記事では扱いません）．

以下では，$\Omega$ と ${\mathcal F}$の具体例と，基本となる実装例を示していきます．

### 見本空間 Ω
**ガチャ**（カプセルおもちゃ販売機）を例にとります（講義でそうしてます）．

ガチャからカプセルを取り出す**試行**を繰り返したところ，次のような種類のマスコットが，ランダムに出てきたとします．$$犬, 鳥, 犬, 人, 人, 鳥, 犬, · · · , 人$$このような**試行結果**の列を，**試行列**といいます．試行結果を自然言語で表すと，紙に書くのが大変なので，適当な記号や数字による，ラベリングがなされます．例えば，次のようにできます．$$1:犬,\quad 2:鳥,\quad 3:人$$個々のラベルを，**見本**といいます．見本を使うと，先の試行列は，次のように簡潔に書けます．$$1, 2, 1, 3, 3, 2, 1, \cdots, 3$$

次に，確率論では，試行結果の全種類を網羅した**見本帳**を作ります．例えば，試行結果が$1,2,3$の3種類であるとき，見本帳を次のような集合で表します．
$$\Omega :=\\{1,2,3\\}$$このような見本帳を表す集合を，**見本空間**といいます．

この$\Omega$は，Sympyの`FiniteSet`を使って，次のように書くことができます．

```python:python3
from sympy import FiniteSet, EmptySet

Om = FiniteSet(1,2,3)
print('Om = ' + str(Om))
```

```console:実行結果
Om = {1, 2, 3}
```

### 事象

続いて，確率論では，見本の様々なクラス分けを導入します．
例えば，哺乳類と2足歩行（おもちゃだけど）に興味があるなら，次のようなクラス$A_1$, $A_2$が導入できます．

| 定義 | 人為的な意味付け |
|:-:|:--:|
| $A_1:=\\{1,3\\}$ | 哺乳類 |
| $A_2:=\\{2,3\\}$ | 2足歩行 |

このようなクラス分けを表す，$\Omega$の部分集合を，**事象**といいます．

同様に，`FiniteSet`を使って次のように書けます．

```python:python3
A1 = FiniteSet(1,3)
A2 = FiniteSet(2,3)
print('A1 = ' + str(A1))
print('A2 = ' + str(A2))
```

```console:実行結果
A1 = {1, 3}
A2 = {2, 3}
```

以上を踏まえて，事象$A_1$, $A_2$を要素とする，**族**（集合を要素とする集合のこと）を考えます．$$F:=\\{A_1,A_2\\}=\\{\\{1,3\\},\\{2,3\\}\\}$$この族は，考察中の事象のリストを与えることになります．

`FiniteSet`の引数は`FiniteSet`でもよいので，次のように書けます．

```python:python3
F = FiniteSet(A1,A2)
print('F = ' + str(F))
```

```console:実行結果
F = {{1, 3}, {2, 3}}
```

### σ加法族

ここからが本題です．
確率論にはルールがありまして，事象の集合演算結果もまた，事象でなければなりません．数学でよく目にするルールですね（ある集合の要素間に演算を導入したとき，演算結果がその集合から出てはいけないという）．

残念ながら，さきほどの事象のリスト$F$は，このルールを満たしていません．例えば，$A_1\cap A_2=\\{3\\}$は，$F$の要素ではありません．
そこで，不完全な$F$を拡張して，完全なる全品リストを得ることを考えます．
そのためのチェックシートが整備されています．

**□定義（σ加法族）**：　$\Omega$の部分集合族$\mathcal F$で，次の3条件を満すものをσ加法族という．

1. $\Omega\in{\mathcal F}$.
2. $A\in {\mathcal F}\implies \overline A\in {\mathcal F}$.
3. $A_1,A_2,\cdots,A_n\in {\mathcal F}\implies A_1\cup A_2\cup\cdots\cup A_n\in {\mathcal F}$.　($n$は任意の自然数)

σ加法族は集合演算で閉じる．□

この定義を満たすように，元の$F$を拡張していけば，完全なる$\mathcal F$，すなわちσ加法族が得られます．
ところが，それを手作業で，間違いなく実行するのは，案外大変なのです．
そこで，次のようなプログラムを書いてみた次第です．

## 実装 ー σ加法族の自動生成

まず，**定義（σ加法族）**を満たすか，チェックする関数を作ります．

```python:python3
def is_sigma_algebra(Om,FF):
    is_ok = True

    # (条件1)のチェック
    is_ok = is_ok and (Om in FF)

    # for文で全要素を巡るためのリスト化
    elm=list(FF) #リスト化
    n=len(elm)   #リストの要素数

    # (条件2)のチェック
    for i in range(0,n):
        is_ok = is_ok and (Om-elm[i] in FF)

    # (条件3)の全数チェック 
    for i in range(0,n):
        for j in range(i,n):
            is_ok = is_ok and (elm[i]+elm[j] in FF)

    return is_ok
```

上述したように，さきほどの`F`は，チェックを通りません．

```python:python3
print( 'F is sigma-algebra: ' + str(is_sigma_algebra(Om,F)))
```

```console:実行結果
F is sigma-algebra: False
```

次に，**定義（σ加法族）**に出てくる集合演算結果（余事象と和事象）を追加する関数を作ります．

```python:python3
#余事象を追加する関数
def append_complements(Om,F):
    result=F

    #いまある全ての事象の余事象を追加する
    for elm in list(F):
        celm = Om - elm
        result = result + FiniteSet(celm)

    return result

#和事象を追加する関数
def append_unions(F):
    result=F
    Flist=list(F)
    Fn=len(Flist)

    #いまある全ての事象の組み合わせの和事象を追加する
    for i in range(0,Fn):
        for j in range(i,Fn):
            result = result + FiniteSet(Flist[i]+Flist[j])

    return result
```

以上を組み合わせて，σ加法族を生成する関数を作ります．

```python:python3
def generate_sigma_algebra(Om,F):
    cur_F = F + FiniteSet(Om)
    prev_F = EmptySet()

    #族が変化しなくなるまでループ
    while(prev_F != cur_F):
        prev_F = cur_F
        #余事象の追加
        cur_F = append_complements(Om,cur_F)
        #和事象の追加
        cur_F = append_unions(cur_F)

    return cur_F
```

さきほどの`F`で生成されるσ代数を求めてみます．

```python:python3
FF = generate_sigma_algebra(Om,F)
print('FF = ' + str(FF))
print( 'FF is sigma-algebra: ' + str(is_sigma_algebra(Om,FF)))
```

```console:実行結果
FF = {EmptySet(), {1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}
FF is sigma-algebra: True
```

$\Omega$のべき集合が出てきました．他の例も試してみます．

```python:python3
A3 = FiniteSet(2)
F2 = FiniteSet(A1,A3)
print('F2 = ' + str(F2))
FF2 = generate_sigma_algebra(Om,F2)
print('FF2 = ' + str(FF2))
print( 'FF2 is sigma-algebra: ' + str(is_sigma_algebra(Om,FF2)))
```

```console:実行結果
F2 = {{2}, {1, 3}}
FF2 = {EmptySet(), {2}, {1, 3}, {1, 2, 3}}
FF2 is sigma-algebra: True
```

今度はもう少し小さいですね．


### 確率変数で生成されるσ加法族

以上の簡単な応用として，確率変数からσ加法族を生成することもできます．下記が実装例です．

+ https://github.com/ktysd/prob-tut/blob/master/Jupyter/Ch02.ipynb

## おわりに

ざざっと実装してみたのですが，素人の自己流なので，もっと良い方法があるような気もします．
例えば，上述のアルゴリズムには２重ループが何度もでてくるので，事象の数が多いと厳しいですよね．
何かお気づきの点があれば，ぜひぜひご教示ください．よろしくお願いいたします．