[Math] 班佛定律 Benford's law

數學知識：班佛定律

簡介：班佛定律的原理

自然界首數的機率並不是平均分配的，像是河川的長度、戶頭的金錢，甚至是請人們隨機選出的四位數字

在 $b$ 進位制中，以數 $n$ 起頭的數出現的機率為
$$\begin{array}{rl}P(n)& =\frac{lo{g}_b(n+1)-lo{g}_b(n)}{lo{g}_b(b)-lo{g}_b(1)}\\ & =lo{g}_b(n+1)-lo{g}_b(n)\\ & =lo{g}_b(\frac{n+1}{n})\\ & =lo{g}_b(1+\frac{1}{n})\end{array}$$

一組平均增長的數據開始時，增長得較慢
由最初的數字 $a$ 增長到另一個數字 $a+1$ 起首的數的時間
必然比 $a+1$ 起首的數增長到 $a+2$，需要更多時間，所以出現率就更高了
從數數目來說，順序從1開始數，1,2,3,...,9，從這點終結的話，所有數起首的機會似乎相同
但 9 之後的兩位數 10 至 19，以 1 起首的數又大大拋離了其他數了
而下一堆 9 起首的數出現之前，必然會經過一堆以 2,3,4,...,8 起首的數
若果這樣數法有個終結點，以 1 起首的數的出現率一般都比 9 大
或者換個角度來看
增長的速度需為幾何成長，才有此現象
因為無論是 $2^x$、${10}^x$、$e^x$，取 log 後，皆會變為線性
而由於每段分配類似梯型，因基數極大，導致上底與下底差異小，所以會由高決定機率
若是 [1,10)，分別機率為 $lo{g}_{10}(2)-lo{g}_{10}(1)$、$lo{g}_{10}(3)-lo{g}_{10}(2)$...$lo{g}_{10}(10)-lo{g}_{10}(9)$
而若是 [10,100)，分別機率為 $lo{g}_{10}(20)-lo{g}_{10}(10)$、$lo{g}_{10}(30)-lo{g}_{10}(20)$...$lo{g}_{10}(100)-lo{g}_{10}(90)$
實際上也是一樣的，因為 $log$ 相減等於相除

十位數首數分配機率

人們隨性所寫出的數，銀行帳簿的存款數據，河流的區域
面積，都市的人口數據，公司的收支數據等，都會有底下的分佈比例

$n$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$
$P(n)$	30.1%	17.6%	12.5%	9.7%	7.9%	6.7%	5.8%	5.1%	4.6%

${10}^x$ 首數分配

import numpy as np
import matplotlib.pyplot as plt
 
def getFirst(x):
    if x < 10:
        return int(x)
    else:
        return getFirst(x//10)
f = np.vectorize(getFirst)
 
data = np.random.rand(10000)
data = np.floor(10**data)
x = f(data)
plt.hist(x, bins=range(1, 11), density=True)
plt.show()

台股股價首數分佈

import requests
import pandas as pd
import time
from pyquery import PyQuery
import matplotlib.pyplot as plt
 
 
# 得到首位數
def getFirst(x):
    if x < 10:
        return int(x)
    else:
        return getFirst(x//10)
 
    
# 得到股票代碼        
def getSymbolList():
    url = "http://isin.twse.com.tw/isin/C_public.jsp?strMode=2"
    dom = PyQuery(url)
    data = []
    for i, row in enumerate(dom('table.h4 tr').items()):
        if i < 2:
            continue
        item = []
        for j, col in enumerate(row.find('td')):
            if col.text:
                if j == 0: # 分開代碼與名字
                    item.extend(col.text.strip().split('\u3000'))
                else:
                    item.append(col.text.strip().replace('\u3000', ''))
        if item:
            data.append(item)
        else: # 只取上市部分
            return data
            
    return data
    
    
# Alpha Vantage API
def alphaVantage(function, symbol, interval='1min'):
     url = "https://www.alphavantage.co/query"
     headers = {"User-Agent": "Mozilla/5.0 (Windows NT 6.1; rv:57.0) Gecko/20100101 Firefox/57.0",
                "Connection": "keep-alive",
               }
     api_key = "xxxx" #your key
 
     params = {"function": function,
             "symbol": symbol,
             "interval": interval,
             "apikey": api_key}
     c = 0
     while True:
         page = requests.get(url, params=params)#, headers=headers)
         data = page.json()
         if len(data) >= 2:
             return data
 
         c += 1
         if c > 5:
             break
         time.sleep(1)
     print(data, page.url)
     return data   
 
     
# 得到股票價格
def getStockPrice(symbol):
    history= alphaVantage('TIME_SERIES_INTRADAY', symbol)
    last = history.get('Meta Data')
    if last is None:
        return 0.0
    last = last['3. Last Refreshed']
    price = history['Time Series (1min)'][last]['4. close']
    if price:
        print(symbol, float(price))
        return float(price)
          
     
if __name__ == '__main__':
    data = getSymbolList()
    df = pd.DataFrame(data)
    df = df.ix[:,[0,1]].dropna()
    df = df.rename(index=int, columns={0:'symbol', 1:'name'})
    df.ix[:, 'symbol'] = df['symbol'] + '.TW'
    df.ix[:, 'price'] = 0
    
    # 取得價格
    df = df.reset_index(drop=True)
    for i, symbol in enumerate(df['symbol']):
        if df.ix[i, 'price'] == 0:
            df.ix[i, 'price'] = getStockPrice(symbol)
        
    df.ix[:, 'first'] = df['price'].apply(getFirst)
    
    plt.hist(df['first'], bins=range(1, 11), density=True)
    plt.show()  

參考

不要相信你的直覺
Wiki 班佛定律
52 班佛法則...不會丟擲骰子的上帝？
A Statistical Derivation of the Signi cant-Digit Law