1～12などの番号が書かれたイス

プレイヤー二人は交互に以下を行う

相手が座るイスを予想し、電流を仕掛ける

電気イスでなければ得点獲得し、イスは撤去される

電流イスに座ると得点没収され、イスはそのまま

最終的にポイントが高い方が勝利

3回電流を喰らうか、40点先取されたら負け

import random

headers = ["ターン数", "選択プレイヤー", "選択数字", "電気プレイヤー", "仕掛数字", "結果", "P1得点", "P2得点", "P1電気回数", "P2電気回数", "残り椅子"]

class ElectricChairGame:
    def __init__(self):
        self.chairs = {i: i for i in range(1, 13)}
        self.player_scores = [0, 0]
        self.electric_fails = [0, 0]
        self.current_player = 0
        self.electrified_chair = None

    def choose_chair(self, chair_number):
        # 選択した数字と電気椅子の数字が一致する場合のみ電気椅子確定
        if chair_number == self.electrified_chair:
            self.electric_fails[self.current_player] += 1
            self.player_scores[self.current_player] = 0
            return "電気椅子"
        else:
            points = self.chairs[chair_number]
            self.player_scores[self.current_player] += points
            del self.chairs[chair_number]
            return "得点獲得"

    def set_electric(self, chair_number):
        self.electrified_chair = chair_number
        return True

    def check_game_end(self):
        message = ["", ""]
        for player in range(2):
            if self.electric_fails[player] >= 3:
                winner = 2 if player == 0 else 1
                message = [winner, "電気椅子3回"]
                return True, message
            if self.player_scores[player] >= 40:
                winner = player + 1
                message = [winner, "スコア40"]
                return True, message

        if len(self.chairs) == 1:
            winner = 1 if self.player_scores[0] > self.player_scores[1] else 2
            message = [winner, "スコア差"]
            return True, message

        return False, message

def simulate_game(echo=False):
    game = ElectricChairGame()
    turn,selecting_player = 1,1
    results = []
    while True:
        # 電気椅子設置
        electric_player = 2 if selecting_player == 1 else 1
        electric_choice = random.choice(list(game.chairs.keys()))
        game.set_electric(electric_choice)

        # 椅子選択
        chair_choice = random.choice(list(game.chairs.keys()))
        game.current_player = selecting_player - 1
        result = game.choose_chair(chair_choice)

        # ログ出力
        log_line = [
            turn,
            selecting_player,
            chair_choice,
            electric_player,
            electric_choice,
            result,
            game.player_scores[0],
            game.player_scores[1],
            game.electric_fails[0],
            game.electric_fails[1],
            sorted(list(game.chairs.keys()))
        ]
        results.append(log_line)
        # ゲーム終了判定
        game_end, message = game.check_game_end()
        if game_end:
          if message[0]==1:
              win_score = game.player_scores[0]
          else:
              win_score = game.player_scores[1]
          message+=[turn ,
            win_score,
            game.player_scores[0],
            game.player_scores[1],
            game.electric_fails[0],
            game.electric_fails[1],
            len(list(game.chairs.keys()))]
          if echo:
            print(", ".join(headers))
            for log_line in results:
              print(log_line)
            print(f"ゲーム終了: 勝利プレイヤー : {message[0]}, 勝因 : {message[1]}")
          return message, results

        # プレイヤー交代
        selecting_player = 2 if selecting_player == 1 else 1
        turn += 1

message, log_line = simulate_game(True)

ターン数, 選択プレイヤー, 選択数字, 電気プレイヤー,
仕掛数字, 結果, P1得点, P2得点, P1電気回数, P2電気回数, 残り椅子

[1, 1, 2, 2, 3, '得点獲得', 2, 0, 0, 0, [1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]]

[2, 2, 11, 1, 5, '得点獲得', 2, 11, 0, 0, [1, 3, 4, 5, 6, 7, 8, 9, 10, 12]]

[3, 1, 6, 2, 7, '得点獲得', 8, 11, 0, 0, [1, 3, 4, 5, 7, 8, 9, 10, 12]]

[4, 2, 1, 1, 8, '得点獲得', 8, 12, 0, 0, [3, 4, 5, 7, 8, 9, 10, 12]]

[5, 1, 3, 2, 10, '得点獲得', 11, 12, 0, 0, [4, 5, 7, 8, 9, 10, 12]]

[6, 2, 12, 1, 7, '得点獲得', 11, 24, 0, 0, [4, 5, 7, 8, 9, 10]]

[7, 1, 8, 2, 4, '得点獲得', 19, 24, 0, 0, [4, 5, 7, 9, 10]]

[8, 2, 9, 1, 7, '得点獲得', 19, 33, 0, 0, [4, 5, 7, 10]]

[9, 1, 10, 2, 7, '得点獲得', 29, 33, 0, 0, [4, 5, 7]]

[10, 2, 7, 1, 5, '得点獲得', 29, 40, 0, 0, [4, 5]]

ゲーム終了: 勝利プレイヤー : 2, 勝因 : スコア40

game_result, game_los = [], []
for i in range(10000):
  message, log_line = simulate_game()
  game_result.append([i+1] + message)
  for log in log_line:
    game_los.append([i+1] + log)

import pandas as pd
cols = ["ゲーム数", "勝利プレイヤー", "勝因","ターン数","獲得点数","P1得点", "P2得点", "P1電気回数", "P2電気回数", "残椅子個数"]
game_result_df = pd.DataFrame(game_result, columns=cols)
game_log_df = pd.DataFrame(game_los, columns=["ゲーム数"] + headers)

game_result_df.head()

h=(d+[26(m+1)/10]+Y+[Y/4]-2[y/100]+[y/400])mod 7

0=土曜日

1=日曜日

2=月曜日

3=火曜日

4=水曜日

5=木曜日

6=金曜日

# モジュロ演算で曜日を求めるコード
def zeller(year, month, day):
    # 1月と2月は前年の13月、14月として扱う
    if month < 3:
        month += 12
        year -= 1

    y = year % 100
    century = year // 100
    
    # ツェラーの公式より
    h = (day + (26 * (month + 1)) // 10 + y + (y // 4) - (2 * century) + (century // 4)) % 7
    
    # 曜日
    # 0=土曜日, 1=日曜日, 2=月曜日, 3=火曜日, 4=水曜日, 5=木曜日, 6=金曜日
    return h

weekdays = ["土曜日", "日曜日", "月曜日", "火曜日", "水曜日", "木曜日", "金曜日"]

# 使用例
year, month, day = 2025, 2, 1

h = zeller(year, month, day)
weekday_name = weekdays[h]
print(f"{year}年{month}月{day}日は{weekday_name}です。")

なつき先輩の誕生日は
1992年7月19日

解説動画はこちら



ベンフォードの法則

今回はベンフォードの法則のお話です

この法則は、自然界に出てくる多くの数値の最初の桁の分布が
この法則によると

math.log10((d+1)/d)

ここから先のコードをGoogle Colabで試す場合は
以下をインストールしてください

pip install japanize_matplotlib

import math
import numpy as np
import scipy.stats as stats
import japanize_matplotlib
import matplotlib.pyplot as plt

# ベンフォードの理論的な確率分布
def benford_distribution():
    return [math.log10(1 + 1 / d) for d in range(1, 10)]

benford_probs = benford_distribution()

# 棒グラフの描画
labels = [1, 2, 3, 4, 5, 6, 7, 8, 9]
bar_width = 0.35
x = np.arange(len(labels))
plt.bar(x - bar_width/2, benford_probs, width=bar_width, label='理論的確率', color='blue')
plt.xlabel('最初の桁')
plt.ylabel('確率')
plt.title('ベンフォードの法則の比較')
plt.xticks(x, labels)
plt.legend()
plt.tight_layout()
plt.show()

ベンフォードの法則は何に使えるのか

次のような不正の調査などで利用されているらしいです

ここから先はどのように不正を発見するのかを
見てみましょう。




・カイ二乗検定を用いる方法

カイ二乗検定は観察されたデータと期待される
データの間の差を評価するための統計的手法のことです
カイ二乗検定を用いたベンフォードの法則比率との比較


ここでは理論値と
観測値として、操作した比率のデータで検証していきます。

以下のコードで描画させる事ができます。

import numpy as np
import scipy.stats as stats
import japanize_matplotlib
import matplotlib.pyplot as plt

# サンプルサイズを設定
n_samples = 1000  # 適宜変更

# ベンフォードの法則に従った最初の桁の理論的確率
benford_probs = [0.301, 0.176, 0.125, 0.097, 0.079, 0.067, 0.058, 0.051, 0.046]

# 観測された出現確率を設定
observed_probs = [0.295, 0.172, 0.166, 0.092, 0.074, 0.062, 0.053, 0.046, 0.040]

# 最初の桁のラベル
labels = [1, 2, 3, 4, 5, 6, 7, 8, 9]

# 棒グラフの幅
bar_width = 0.35
x = np.arange(len(labels))

# 棒グラフの描画
plt.bar(x - bar_width/2, benford_probs, width=bar_width, label='理論的確率', color='blue')
plt.bar(x + bar_width/2, observed_probs, width=bar_width, label='観測された確率', color='orange')

# グラフの設定
plt.xlabel('最初の桁')
plt.ylabel('確率')
plt.title('ベンフォードの法則の比較')
plt.xticks(x, labels)
plt.legend()
plt.tight_layout()
plt.show()

数字の3が出てくる確率を多くしました。
これが理論上の比率と合うのかを検定します。

# 出現数を計算
observed_counts = [int(p * n_samples) for p in observed_probs]
expected_counts = [int(p * n_samples) for p in benford_probs]

print("観測数 : ",observed_counts)
print("期待値 : ",expected_counts)

# カイ二乗検定を実施
chi2_stat, p_value = stats.chisquare(f_obs=observed_counts, f_exp=expected_counts)

# 結果を表示
print("カイ二乗統計量:", chi2_stat)
print("p値:", p_value)

# 結論
alpha = 0.05  # 有意水準
if p_value < alpha:
    print("ベンフォードの法則に従わないと判断される")
else:
    print("ベンフォードの法則に従うと判断される")

この操作された比率の場合は
ベンフォードの法則に従わないと判断されるようです。





理論上の比率と観測された比率が正しいのか
異なるのかを仮説に置き
一般的には帰無仮説 には差がない(比率と等しい)とします。



となり、閾値より大きくなった場合は
P値(確率)が低くなり

あり得ない事が起きているとして帰無仮説を棄却し
対立仮説 : 観察された比率は期待される比率と異なる
を採択することになります。


カイ二乗分布の検定統計量の表を作ってみましょう

import numpy as np
import pandas as pd
from scipy.stats import chi2

# 有意水準
alpha_levels = [0.99, 0.975, 0.95, 0.9, 0.1, 0.05, 0.025, 0.01]
# 自由度
degrees_of_freedom = range(1, 10)

# カイ二乗分布の臨界値を計算
chi_square_table = {'自由度/有意水準': degrees_of_freedom}

for alpha in alpha_levels:
    chi_square_table[f'α = {alpha}'] = [chi2.ppf(1 - alpha, df) for df in degrees_of_freedom]

# カイ二乗値のDataFrameを作成
chi_square_df = pd.DataFrame(chi_square_table)
chi_square_df

こんな感じで、有意水準と自由度で
閾値を求められます。

紙で行う場合は、この表を見て行いますが
Pythonプログラミングでは、これを使わず
直接確率を求められます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2

# 自由度
df = 8
# 有意水準
alpha = 0.05
# カイ二乗の臨界値を計算
critical_value = chi2.ppf(1 - alpha, df)

# 観測したカイ二乗値
observed_chi_squared = sum([(o*1000-b*1000)**2/(b*1000) for b,o in zip(benford_probs,observed_probs)])

# xの範囲を設定
x = np.linspace(0, 30, 1000)
# カイ二乗分布の確率密度関数を計算
y = chi2.pdf(x, df)

# グラフを描画
plt.figure(figsize=(10, 6))
plt.plot(x, y, label=f'Chi-squared Distribution (df={df})', color='blue')
plt.fill_between(x, y, where=(x >= critical_value), color='red', alpha=0.5, label='Rejection Region (p < 0.05)')
plt.axvline(critical_value, color='red', linestyle='--', label=f'Critical Value = {critical_value:.2f}')
plt.axvline(observed_chi_squared, color='green', linestyle='--', label=f'Observed Chi-squared = {observed_chi_squared:.2f}')
plt.title('Chi-squared Distribution with 5% Significance Level')
plt.xlabel('Chi-squared Value')
plt.ylabel('Probability Density')
plt.legend()
plt.grid()
plt.show()

この場合は求めたカイ二乗値は
閾値を超える値になるため
P値は低くなり、対立仮説の採択、となります。





データのリンク先

こちらをColabで行う場合は
ファイル置き場においてください。


データを読み込むコードはこれです。

import pandas as pd
import numpy as np

file_path = "major_results_2020.xlsx"
  # 8行目をヘッダーとして指定（0-indexed）
df = pd.read_excel(file_path, header=8)

df.iloc[0:2, 35:]

これを集計して最初の数値のカウントと
比率を求めます。

data = df.iloc[:, 35:]

# 先頭の数字をカウントするための関数
def count_leading_digits(series):
    leading_digits = series.astype(str).str[0]
    return leading_digits.value_counts()

# 全列のデータを結合して先頭の数字をカウント
all_leading_digits = data.astype(str).stack().str[0]
digit_counts = all_leading_digits.value_counts().sort_index()
digit_counts = digit_counts[~digit_counts.index.str.contains('-')]

# 結果を表示
print(digit_counts)

# 比率を計算
total_count = digit_counts.sum()
leading_digit_ratios = digit_counts / total_count
print(leading_digit_ratios)

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.stats import norm
from collections import Counter

# ベンフォードの法則による理論的な確率分布
def benford_distribution():
    return [math.log10(1 + 1 / d) for d in range(1, 10)]

# ベンフォードの理論的分布
benford_dist = benford_distribution()

# サンプルデータの分布
sample_distribution = list(leading_digit_ratios.values)

# グラフ描画
x = range(1, 10)
fig, ax = plt.subplots(figsize=(10, 6))
ax.bar(x, sample_distribution, width=0.4, label="Sample Distribution", align="center", alpha=0.7)
ax.bar([xi + 0.4 for xi in x], benford_dist, width=0.4, label="Benford Distribution", color="orange", alpha=0.7)

# グラフの設定
ax.set_xticks(x)
ax.set_xticklabels([str(d) for d in x])
ax.set_xlabel("Leading Digit")
ax.set_ylabel("Frequency")
ax.set_title("Comparison of Sample Distribution and Benford's Law")
ax.legend()
plt.show()

カイ二乗検定を行うと

# カイ二乗検定の実施
import scipy.stats as stats

# サンプルの頻度を期待値に基づいて計算
expected_frequencies = [sum(sample_distribution) * p for p in benford_dist]
observed_frequencies = sample_distribution

# カイ二乗検定を実行
chi2_stat, p_value = stats.chisquare(observed_frequencies, f_exp=expected_frequencies)

# 結果を表示
print("カイ二乗統計量:", chi2_stat)
print("p値:", p_value)

# 有意水準を設定
alpha = 0.05
if p_value < alpha:
    print("帰無仮説を棄却します。サンプルデータはベンフォードの法則に従っていない可能性があります。")
else:
    print("帰無仮説を棄却できません。サンプルデータはベンフォードの法則に従っている可能性があります。")

このデータの場合は
ベンフォードの法則に従っていないとは言えないようですね

正しいデータである可能性が高いです。



まとめ

会計データなどがある場合は
これで比率が正しいかを試してみると
面白いかもしれませんね。

今日はこれまでです
それでは