Zenn
🦋

あみだくじの問題をjuliaで！（その2）

2024/08/07に公開
Julia
math
確率
あみだくじ
idea

はじめに

Xのポストで@yusui_mathさんがこのような数学の問題が投稿されました。

https://x.com/yusui_math/status/1817484461841515005

@yusui_mathさんはMathlogのサイトで表現論・群論を用いて解いている様子が公開されています。

https://mathlog.info/articles/xknnpdUZFn8ww5Hkufov

以前，私もZennのサイトであみだくじの問題を考えたことがあったので，今回もJuliaで考えてみることにしました。

https://zenn.dev/dannchu/articles/c8835c0eb8caa2

まずは実験です

  • あみだくじの横棒をn回を引く。（デフォルトはn=10^4）
  • 元の位置に戻る（恒等置換）かチェックする。
  • これをm回行って割合を調べる。（デフォルトはm=10^4）
  • この実験を10回繰り返す。
  • パッケージは抽象代数パッケージAbstractAlgebra.jlを使いました。
#抽象代数パッケージ
using AbstractAlgebra 

#3つの置換
σ₁ = perm([2,1,3,4])
σ₂ = perm([1,3,2,4])
σ₃ = perm([1,2,4,3])
X = [σ₁ , σ₂ , σ₃ ]

# n回横線を引いて，元の位置に戻るかチェックする。これをm回行って割合を調べる。
function prob(n=10^4,m=10^4)
    k = 0
    for i = 1:m
        f =  perm([1,2,3,4])
        for i = 1:n
            f *= rand(X)
        end
        if f == perm([1,2,3,4])
            k += 1
        end
    end
    k/m
end

#割合を10回調べてみる。1/12 = 0.08333 に近ければいい感じ
for _=1:10
    @show prob()
end

prob() = 0.0872
prob() = 0.083
prob() = 0.0835
prob() = 0.079
prob() = 0.0882
prob() = 0.0898
prob() = 0.0793
prob() = 0.0845
prob() = 0.0812
prob() = 0.0831

遷移行列を求めて，確率を計算する。

  • あみだくじの横線の数nによって，求める確率がどのようになるかを調べます。
  • 元に戻るケースだけでなく，全てのケースの確率を考えます。
  • 全ての確率のn-1とnの関係を表す遷移行列Aを求めます。

S(4)の置換群

今回準備した，3つの置換σ₁ = perm([2,1,3,4])，σ₂ = perm([1,3,2,4])，σ₃ = perm([1,2,4,3])の合成で，4!=24通りの置換は全て表すことができるようです。6回の合成で全てで揃いました。

σ₁ = Perm([2,1,3,4]);
σ₂ = Perm([1,3,2,4]);
σ₃ = Perm([1,2,4,3]);
X = [A,B,C]

for j =1:5
    for i = 1:length(X)
        push!(X,X[i]*σ₁,X[i]*σ₂,X[i]*σ₃)
    end

    X = X |> unique
    println("合成$(j+1)回で置換は$(X|>length)種類")
end

合成2回で置換は9種類
合成3回で置換は15種類
合成4回で置換は20種類
合成5回で置換は23種類
合成6回で置換は24種類

遷移行列Aを求めます。

あみだくじにn本の横線を入れた時に，置換\sigma_k\quad(k=1,2,...,24)となる確率を\sigma_k(n)とします。

\begin{pmatrix} \sigma_1(n)\\ \sigma_2(n)\\ \vdots\\ \sigma_{24}(n)\\ \end{pmatrix} = A \begin{pmatrix} \sigma_1(n-1)\\ \sigma_2(n-1)\\ \vdots\\ \sigma_{24}(n-1)\\ \end{pmatrix}

上の式を満たす24\times 24の正方行列Aを求めます。

Juliaのコード

3つの置換と，24通りのS(4)を準備

  • 24通りの置換を作成するのに，組み合わせパッケージCombinatorics.jlを利用しました。
  • 1番目の\sigma_1が恒等置換です。
using AbstractAlgebra,Combinatorics

# 定義された置換
σ₁ = Perm([2,1,3,4]);
σ₂ = Perm([1,3,2,4]);
σ₃ = Perm([1,2,4,3]);

# 置換群S4のすべての置換
X = Perm.(permutations([1,2,3,4]))

24-element Vector{Perm{Int64}}:
()
(3,4)
(2,3)
(2,3,4)
(2,4,3)
(2,4)
(1,2)
(1,2)(3,4)
(1,2,3)
(1,2,3,4)
⋮
(1,3,4)
(1,3)(2,4)
(1,3,2,4)
(1,4,3,2)
(1,4,2)
(1,4,3)
(1,4)
(1,4,2,3)
(1,4)(2,3)

遷移行列を求めるコード

A，B，Cをかけて得られた置換がXの何番目に当たるかを調べる方法がちょっとわからなかったので，ChatGPT-4oにも助けてもらいました。 findfirst(p -> p == next_perm, X)と書くようですね。


# 各置換に対して行列を生成
function generate_transition_matrix(X, A, B, C)
    n = length(X)
    transition_matrix = zeros(Rational{Int}, n, n)
    
    for i in 1:n
        current_perm = X[i]
        next_perms = [current_perm * A, current_perm * B, current_perm * C]
        
        for next_perm in next_perms
            j = findfirst(p -> p == next_perm, X)
            if j !== nothing
                transition_matrix[i, j] += 1 // 3
            end
        end
    end
    return transition_matrix
end


# 置換のリストに基づいて遷移行列を生成
transition_matrix = generate_transition_matrix(X, σ₁ , σ₂ , σ₃) 

# 結果を表示
println("Transition Matrix:")
for row in eachrow(transition_matrix)
    println(row)
end

Transition Matrix:
Rational{Int64}[0, 1//3, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[1//3, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[1//3, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[0, 0, 1//3, 0, 0, 1//3, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[0, 0, 0, 1//3, 1//3, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[1//3, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[0, 1//3, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0]
Rational{Int64}[0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0]
Rational{Int64}[0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 1//3, 0]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 1//3]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 0, 0, 1//3, 1//3, 0, 0, 0]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 0, 1//3, 0, 0, 1//3, 0, 0]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 0, 1//3]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 0, 0, 1//3]
Rational{Int64}[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1//3, 0, 0, 0, 1//3, 1//3, 0]

値をチェック！

  • n=1,2,...,10で，確率を調べてみました。求めるのは\sigma_1(n)です。
\begin{pmatrix} \sigma_1(n)\\ \sigma_2(n)\\ \vdots\\ \sigma_{24}(n)\\ \end{pmatrix} = A^n \begin{pmatrix} 1\\ 0\\ \vdots\\ 0\\ \end{pmatrix}
K = zeros(Rational{Int},24)
K[1] = 1
for i = 1:10
    println("$(i)回で元に戻る確率は$((transition_matrix^i * K)[1])")
end

1回で元に戻る確率は0//1
2回で元に戻る確率は1//3
3回で元に戻る確率は0//1
4回で元に戻る確率は17//81
5回で元に戻る確率は0//1
6回で元に戻る確率は115//729
7回で元に戻る確率は0//1
8回で元に戻る確率は283//2187
9回で元に戻る確率は0//1
10回で元に戻る確率は6643//59049

@yusui_mathさんのページの確率と一致しました。

nが偶数の場合の極限は\frac1{12}のようなので，\sigma_1(30)を調べてみました。

(transition_matrix^30 * K)[1] |> Float64 , 1/12

(0.08370280962869549, 0.08333333333333333)

いい感じです！

Discussion