解説動画はこちら



ベンフォードの法則

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

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

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("帰無仮説を棄却できません。サンプルデータはベンフォードの法則に従っている可能性があります。")

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

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



まとめ

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

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

問題

5リットル入るバケツと3リットル入るバケツがある

この2つだけを使って4リットルの水を測りたい

どうすれば良いでしょうか

ヒント : 水は捨てたり移動したりして良い

幅優先探索(breadth-first search, BFS)を
用いて手順の探索を行っていきます。

バケツをA, B、測りたい水量をXとし
A,Bの水量を記録していきます

途中、以下の6つの状態でうまく切り替えて手順を記録し
どこかが水量Xになったら終了になります。

1. Aを満たす

2. Bを満たす

3. Aを空にする

4. Bを空にする

5. AからBに移す

6. BからAに移す

from collections import deque
from math import gcd

# 幅優先探索(breadth-first search, BFS)を用いて手順の探索を行う
def min_steps_to_measure(A, B, X):
    # 目標の水量がAとBの最大公約数で割り切れない場合は終了
    if X > max(A, B) or X % gcd(A, B) != 0:
        return []

    # キューとSETと状態の親を初期化
    queue, visited, parent = deque(), set(), {}

    # 初期状態をキューに追加
    queue.append((0, 0))  # (容器Aの水量, 容器Bの水量)
    parent[(0, 0)] = None  # 初期状態の親はNone

    while queue:
        a, b = queue.popleft()

        # 目標の水量に達した場合
        if a == X or b == X:
            # 結果を表示
            steps = []
            while (a, b) is not None:
                steps.append((a, b))
                # 親が存在する場合のみアンパック
                if parent[(a, b)] is not None:
                    a, b = parent[(a, b)]
                else:
                    break
            steps.reverse()
            return steps

        # 訪れた状態を記録
        if (a, b) in visited:
            continue
        visited.add((a, b))

        # 次の状態をキューに追加
        # 1. 容器Aを満たす
        if (A, b) not in visited:
            queue.append((A, b))
            parent[(A, b)] = (a, b)
        # 2. 容器Bを満たす
        if (a, B) not in visited:
            queue.append((a, B))
            parent[(a, B)] = (a, b)
        # 3. 容器Aを空にする
        if (0, b) not in visited:
            queue.append((0, b))
            parent[(0, b)] = (a, b)
        # 4. 容器Bを空にする
        if (a, 0) not in visited:
            queue.append((a, 0))
            parent[(a, 0)] = (a, b)
        # 5. 容器Aから容器Bに移す
        transfer_to_B = min(a, B - b)
        if (a - transfer_to_B, b + transfer_to_B) not in visited:
            queue.append((a - transfer_to_B, b + transfer_to_B))
            parent[(a - transfer_to_B, b + transfer_to_B)] = (a, b)
        # 6. 容器Bから容器Aに移す
        transfer_to_A = min(b, A - a)
        if (a + transfer_to_A, b - transfer_to_A) not in visited:
            queue.append((a + transfer_to_A, b - transfer_to_A))
            parent[(a + transfer_to_A, b - transfer_to_A)] = (a, b)

    return []  # 目標の水量を測ることができない場合

# 容器の容量と目標の水量を指定
A = 5  # 容器Aの容量
B = 3  # 容器Bの容量
X = 4  # 測りたい水量

# 最小手数を計算
steps = min_steps_to_measure(A, B, X)

# 結果を表示
if steps:
    print("途中の過程:")
    for step in steps:
        print(f"容器A: {step[0]}L, 容器B: {step[1]}L")
else:
    print("目標の水量を測ることができませんでした。")

途中の過程:

容器A: 0L, 容器B: 0L

容器A: 5L, 容器B: 0L

容器A: 2L, 容器B: 3L

容器A: 2L, 容器B: 0L

容器A: 0L, 容器B: 2L

容器A: 5L, 容器B: 2L

容器A: 4L, 容器B: 3L



このように探索の結果が表示されます。

文字だけだと分かりづらいので
描画してみることにしましょう。



手順を描画するコード

ipywidgetsでスライダーを作り
stepを可変で描画できます。

A,B,Xを色々な設定に変えて
色々遊んでみてください。

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
import ipywidgets as widgets
from IPython.display import display
%matplotlib inline

A = 5  # 容器Aの容量
B = 3  # 容器Bの容量
X = 4  # 測りたい水量

steps = min_steps_to_measure(A, B, X)

# 描画関数
def draw_graph(step):
    a_value = steps[step][0]
    b_value = steps[step][1]
    a_frame_height = A - a_value
    b_frame_height = B - b_value

    # 描画
    fig, ax = plt.subplots(figsize=(8, 4))
    ax.bar(0, a_value, width=0.4, label='Container A Volume (L)', color='blue', linewidth=2, edgecolor='black')
    rect_a = patches.Rectangle((-0.2, a_value), 0.4, a_frame_height, linewidth=2, edgecolor='black', facecolor='none')
    ax.add_patch(rect_a)
    ax.text(0, max(A,B)+1, str(a_value), ha='center', va='bottom', fontsize=12, color='green')
    ax.bar(1, b_value, width=0.4, label='Container B Volume (L)', color='red', linewidth=2, edgecolor='black')
    rect_b = patches.Rectangle((0.8, b_value), 0.4, b_frame_height, linewidth=2, edgecolor='black', facecolor='none')
    ax.add_patch(rect_b)
    ax.text(1, max(A,B)+1, str(b_value), ha='center', va='bottom', fontsize=12, color='green')

    # タイトルとラベルの設定
    ax.set_title(f'Container A : {A} and B: {B} Target Volume : {X} with Frames step : {step}')
    ax.set_xlabel('Containers')
    ax.set_ylabel('Volume (L)')
    ax.set_xticks([0, 1])
    ax.set_xticklabels(['Container A', 'Container B'])
    ax.set_ylim(0, max(A,B)+3)
    plt.show()

# スライダーの作成
step_slider = widgets.IntSlider(value=0, min=0, max=len(steps)-1, step=1, description='Step:')
interactive_plot = widgets.interactive(draw_graph, step=step_slider)
display(interactive_plot)

今回は2つのバケツを使って水量を測るやつを
プログラムにしてみました。

動画ではもう少し手順のかかるやつもやっているので
興味がある場合は見てみてください。

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

年収500万円の人は、日本の年収ランキングで
何位くらいになるでしょうか?

参考 : 日本の人口 1.245億 (2023年)

今回は日本人の
人口が12450万人

平均年収 450万円

年収中央値を 350万円
で設定して行います。

対数正規分布と
おおよその順位を計算するコードです

Google Colabなどに貼り付けて
incomeのスライドを変えてみてください。

import numpy as np
import matplotlib.pyplot as plt
import math
import ipywidgets as widgets
from IPython.display import display
%matplotlib inline

# 対数正規分布のパラメータ設定
median_income = 350  # 中央値（単位: 万円）
mean_income = 450    # 平均値（単位: 万円）
sigma = math.sqrt(2 * math.log(mean_income / median_income)) # 標準偏差を計算
mu = np.log(median_income) # 対数正規分布の μ を計算

# x軸の範囲（年収100万円 ~ 3000万円）
x = np.linspace(100, 3000, 1000)

# パーセンタイルランクを計算する関数
def calc_percentile_rank(income):
    return 1 - 0.5 * (1 + np.sign(income - median_income) * np.sqrt(1 - np.exp(-((np.log(income) - mu) ** 2) / (2 * sigma ** 2))))

# プロットを更新する関数
def update_plot(income):
    # 確率密度関数 (PDF) を計算
    pdf = (1 / (x * sigma * np.sqrt(2 * np.pi))) * np.exp(-((np.log(x) - mu) ** 2) / (2 * sigma ** 2))
    percentile_rank = calc_percentile_rank(income) # ランク計算
    total_population = 124_500_000 # 日本の人口
    rank = (percentile_rank) * total_population
    plt.clf()
    
    # 確率密度関数をプロット
    plt.figure(figsize=(12, 4))
    plt.plot(x, pdf, label="Log-normal Distribution", color="blue")
    plt.axvline(income, color="red", linestyle="--", label=f"Income = {income}") # incomeの位置に線を引く
    plt.axvline(median_income, color='orange', linestyle='--', label=f"Median: {median_income}")
    plt.axvline(mean_income, color='green', linestyle='--', label=f"Mean: {mean_income}")
    plt.fill_between(x, pdf, where=(x >= income), color="lightblue", alpha=0.5) # income以上の範囲を薄い青で塗りつぶす
    plt.title(f"nensyu {income} : Top {percentile_rank*100:.4f}% , Rank {rank:.0f}", fontsize=14)
    plt.xlabel("Income", fontsize=12)
    plt.ylabel("Density", fontsize=12)
    plt.legend()
    plt.grid(alpha=0.5)
    plt.show()

# スライダーを作成
income_slider = widgets.IntSlider(value=500, min=100, max=3000, step=1, description='Income:')
interactive_plot = widgets.interactive(update_plot, income=income_slider)
display(interactive_plot)