Benford's Law Explained — Complete Guide: The Surprising Mathematical Pattern in Real Data

import math
from collections import Counter
from scipy import stats
import pandas as pd

def benfords_expected(digit: int) -> float:
    """Expected frequency of leading digit per Benford's Law"""
    return math.log10(1 + 1/digit)

def get_leading_digit(n) -> int | None:
    """Extract the first significant (non-zero) digit"""
    s = str(abs(float(n))).lstrip('0').replace('.', '')
    for c in s:
        if c.isdigit() and c != '0':
            return int(c)
    return None

def analyze_benford(numbers: list, name: str = "Dataset") -> dict:
    """
    Analyze a numeric dataset against Benford's Law.
    Returns: dict with digit frequencies, deviations, and chi-squared test result.
    """
    digits = [get_leading_digit(n) for n in numbers]
    digits = [d for d in digits if d is not None]
    total = len(digits)

    if total < 300:
        print(f"⚠️  Warning: only {total} values — need 1000+ for reliable Benford analysis")

    counts = Counter(digits)
    results = []

    print(f"\nBenford's Law Analysis: {name}")
    print(f"Sample size: {total:,} numbers")
    print(f"{'Digit':>5} {'Observed':>10} {'Expected':>10} {'Deviation':>12} {'Flag':>15}")
    print("-" * 60)

    observed_freqs = []
    expected_freqs = []

    for d in range(1, 10):
        observed = counts.get(d, 0) / total
        expected = benfords_expected(d)
        deviation = (observed - expected) / expected * 100

        observed_freqs.append(counts.get(d, 0))
        expected_freqs.append(expected * total)

        flag = ""
        if abs(deviation) > 25:
            flag = "⚠️ SUSPICIOUS"
        elif abs(deviation) > 15:
            flag = "⚡ ELEVATED"

        print(f"{d:>5} {observed:>10.1%} {expected:>10.1%} {deviation:>+10.1f}%  {flag}")
        results.append({'digit': d, 'observed': observed, 'expected': expected, 'deviation': deviation})

    # Chi-squared test for overall conformance
    chi2, p_value = stats.chisquare(observed_freqs, expected_freqs)
    conforming = p_value > 0.05

    print(f"\nChi-squared: {chi2:.2f}, p-value: {p_value:.4f}")
    print(f"Benford conformance: {'✅ PASS' if conforming else '❌ FAIL (p < 0.05)'}")

    return {'results': results, 'chi2': chi2, 'p_value': p_value, 'conforming': conforming}

# Example 1: Natural data — company revenues (should conform)
import random
revenues = [random.lognormal(8, 2) for _ in range(5000)]  # log-normal = Benford-conforming
analyze_benford(revenues, "Company Revenues (Natural)")

# Example 2: Fabricated expense data (humans avoid digit 1 when making up numbers)
fabricated = [random.choice([3000, 3200, 3500, 4100, 5200, 6300, 7100, 8200]) +
              random.randint(-100, 100) for _ in range(500)]
analyze_benford(fabricated, "Expense Claims (Suspicious — too many 3-8 starts)")