---

1.
你算錯了，你把你的每一個vlist[n]印出來，
加起來不等於1。程式裡有bug。

我試着把每一項列出來：
(0.210526-0.301030)**2/0.301030 = 0.027210
(0.315789-0.176091)**2/0.176091 = 0.110827
(0.052632-0.124939)**2/0.124939 = 0.041847
(0.105263-0.096910)**2/0.096910 = 0.000720
(0.105263-0.079181)**2/0.079181 = 0.008591
(0.105263-0.066947)**2/0.066947 = 0.021930
(0.052632-0.057992)**2/0.057992 = 0.000495
(0.000000-0.051153)**2/0.051153 = 0.051153
(0.052632-0.045757)**2/0.045757 = 0.001033

chi2(k)加起來乘於19 等於 5.012
p-value = 0.756
在合理範圍。

2.
那個公式是給N很大的時候用的。
chi2 = N * (p-b)^2 / b，
只有N很大的時候分母才可以用 b。

這個公式其實就是（省略sum符號）：
chi2 = (N*p-N*b)^2 / 誤差^2。
當N很大時，誤差=sqrt(N*b)。

可是在本例裡面，9出現的機率的期望值是0.045757
樣本19，因此9出現的次數的期望值是0.87次，不到一次。
這時候只能用Poisson誤差，不能用常態分佈的誤差。

由於Poisson誤差比常態分佈大，真正的chi2值會比上面的

5.012要小，真正的p-value會比0.756要大。

3.
就跟hugoyu一樣，你使用的一系列數字不是
獨立變量，而是下一個數字是前一個數字的累加。
我懷疑不是獨立變量可以用這個定律。

你應該用hugoyu提共的數字的第二列，那才是獨立變量。
用那一組數字得到的
chi2是6.68，p-value是0.57，
依然在合理的範圍內。



※ 引述《Glamsight (安穩殘憶)》之銘言：
: 生活中總是會遇到許許多多數據
: 也不知道究竟是真是假
: 所以現在檢定的時間到了
: 對於這種增量正比於存量的數據通常我們會使用一種神祕的定律 Benford's law
: 簡單來說，該定律描述數據其第一個數據之分布情形
: 具體來說 (10 進位下)，認為數據第一位數字之機率有如下表
:         第一位數字 | 出現機率
:         ------------------------
:         1          | 0.301029996
:         2          | 0.176091259
:         3          | 0.124938737
:         4          | 0.096910013
:         5          | 0.079181246
:         6          | 0.06694679
:         7          | 0.057991947
:         8          | 0.051152522
:         9          | 0.045757491
: 如，各國的地區人口數、各國 GDP、國土面積、放射性元素之半衰期等
: 其數據之第一位數字都成這樣的機率分布，十分之神奇
: 李永樂老師有視頻介紹，並附有該情況下的證明
: https://www.youtube.com/watch?v=CCo4k9Ax7cM
: 本文就考量下一次感染的人數亦可與現有已感染人數呈正比關係
: 及根據 hugoyo 提供之數據嘗試與 Benford's law 進行比較與檢定
: 首先，我們畫張圖
: https://s1.imgs.cc/img/aaaabbvbw.png?_w=750
: 看起來有點像，又有點不像
: 所以我們考慮具體一點的比較，即統計檢定
: 我們善意的假設數據沒有造假，故如下設計
: 虛無假設 H0: 數據沒有造假
:              即，該數據應與 Benford's law 分布相近
: 對立假設 H1: 數據存在造假
:              即，該數據與 Benford's law 分布不相似
: 檢驗方式邏輯大概如下
: 在善意考量沒有造假下，官方數據應該與 Benford's law 相似
: 那麼我們詢問在 Benford's law 的情況下，出現像是這次官方數據的機率大概有多高
: 具體的算法是使用 Chi-Square Test
: 方法來自 Google 文獻
: https://evoq-eval.siam.org/Portals/0/Publications/SIURO/Vol1_Issue1/Testing
: _for_the_Benford_Property.pdf
: (有斷行，請注意)
: 將本次事件之數據第一個數字 k 之機率定為 p(k)
: 該 k 之數字於 Benford's law 下發生機率定為 b(k)
: 則根據文獻可以知道一 Chi-Square Test 檢定量如下
:         檢定量 = 資料筆數 * sum[(p(k)-b(k))^2/b(k)]      (k=1~9)
: Chi-Square Test 自由度為 9-1=8
: 可以知道檢定量值約為 38.35
: 該檢定 p-value = 6.47 * 10^-6 約等於 0
: 也就是說不論設定多嚴苛的顯著水準
: 都會取得拒絕虛無假設的結論 (reject H0)
: 也就是說，我們可以理性推論該數據存在疑慮。
: 以上如果有錯麻煩跟我說。
: 謝謝
: 下附 Python 程式碼 30 行
: def getHeader(value,base,onoff=True):
:     v=int(value/base)
:     if v!=0 and onoff:
:         v=getHeader(v,base,onoff=True)
:     else:
:         v=value
:     return v
: from math import log
: def benfordDistribution(n,base):
:     pr=log(n+1,base)-log(n,base)
:     return pr
: from scipy.stats import chi2
: def getChi(vlist,base):
:     vsum=0
:     for n in range(base-1):
:         pr=benfordDistribution(n+1,base)
:         vsum=vsum+((vlist[n]-pr)**2)/pr
:     v=len(vlist)*vsum
:     pvalue=1-chi2.cdf(v,base-2)
:     return {'Chi-Square 檢定量':v,'p-value':pvalue}
: base=10
: vlist=[45,62,201,218,320,543,639,1356,2033,2744,4515,5974,7711,9692,11794,
: 14384,17213,20448,24335]
: v1list=[getHeader(v,base)/len(vlist) for v in vlist]
: v=getChi(v1list,base)
: for key in v:
:     print(key,v[key])
: ※ 引述《hugoyo (阿佑)》之銘言：
: : 大概大概啦，他們公布的東西就是這樣
: : 天數   日期     確診    累積確診    死亡  累積死亡   死亡率
: : 1     1月18日     45        45       2        2       4.44%
: : 2     1月19日     17        62       0        2       3.23%
: : 3     1月20日    139       201       1        3       1.49%
: : 4     1月21日     17       218       0        3       1.38%
: : 5     1月22日    102       320       3        6       1.88%
: : 6     1月23日    223       543      11       17       3.13%
: : 7     1月24日     96       639       0       17       2.66%
: : 8     1月25日    717      1356      24       41       3.02%
: : 9     1月26日    677      2033      15       56       2.75%
: : 10    1月27日    711      2744      24       80       2.92%
: : 11    1月28日   1771      4515      26      106       2.35%
: : 12    1月29日   1459      5974      26      132       2.21%
: : 13    1月30日   1737      7711      38      170       2.20%
: : 14    1月31日   1981      9692      43      213       2.20%
: : 15    2月 1日   2102     11794      46      259       2.20%
: : 16    2月 2日   2590     14384      45      304       2.11%
: : 17    2月 3日   2829     17213      57      361       2.10%
: : 18    2月 4日   3235     20448      64      425       2.08%
: : 19    2月 5日   3887     24335      65      490       2.01%
: : 猜猜看明天是不是 2.00% 左右??
: : 從一月底開始都穩穩地控制了，死亡率漸漸變低，這樣大家才會比較安心
: : 那麼，明天會多少確診呢?  既然大家覺得數據都是Garbage，再怎麼分析都是Garbage out
: : 還是可以猜猜看明天的Garbage長什麼樣子呀~
: : 不想太認真用分佈的模型來積分成error function
: : 用大家都能懂得多項式就好，抓個3次方。
: : https://i.imgur.com/pqua02g.png  還蠻平的~
: : https://i.imgur.com/ni4FgcQ.png
: : 所以明天中國大概 確診 28297 人，增加 4046 人
: :                  死亡   566 人，增加   76 人
: : 個人希望他們不要真的這樣玩數據。

--
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 203.77.33.103 (臺灣)
※ 文章代碼(AID): #1UF5gW0O (Gossiping)
※ 文章網址: https://www.ptt.cc/bbs/Gossiping/M.1581013664.A.018.html