※ 本文為 zbali.bbs. 轉寄自 ptt.cc 更新時間: 2020-02-06 23:26:49
看板 Gossiping
作者 標題 Re: [問卦] 武漢肺炎也會數學?
時間 Thu Feb 6 23:18:41 2020
生活中總是會遇到許許多多數據
也不知道究竟是真是假
所以現在檢定的時間到了
對於這種增量正比於存量的數據通常我們會使用一種神祕的定律 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
![](https://img.youtube.com/vi/CCo4k9Ax7cM/0.jpg)
本文就考量下一次感染的人數亦可與現有已感染人數呈正比關係
及根據 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 (阿佑)》之銘言:
: ※ 引述《terrymoon (說好的幸福呢)》之銘言:
: : 剛剛在FB看到的
: : https://i.imgur.com/bNLvEbB.jpg
![[圖]](https://i.imgur.com/bNLvEbB.jpg)
: : 170死/7821確診=2.1%
: : 1/31
: : 213死/9800確診=2.1%
: : 2/1
: : 259死/11880確診=2.1%
: : 2/2
: : 304死/14401確診=2.1%
: : 2/3
: : 361死/17238確診=2.1%
: : 怎麼死亡率都剛好在2.1%上下徘徊?
: 大概大概啦,他們公布的東西就是這樣
: 天數 日期 確診 累積確診 死亡 累積死亡 死亡率
: 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/pqua02gh.png)
![[圖]](https://i.imgur.com/ni4FgcQ.png)
: 死亡 566 人,增加 76 人
: 個人希望他們不要真的這樣玩數據。
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.112.90.228 (臺灣)
※ 文章代碼(AID): #1UF2vK-E (Gossiping)
※ 文章網址: https://www.ptt.cc/bbs/Gossiping/M.1581002324.A.F8E.html
推 : 恩恩 跟我想的一樣1F 180.176.0.145 台灣 02/06 23:19
※ 編輯: Glamsight (140.112.90.228 臺灣), 02/06/2020 23:19:35推 : 對!我也這麼覺得!!2F 36.232.18.241 台灣 02/06 23:19
推 : 這個求證有厲害到推一個3F 140.113.170.68 台灣 02/06 23:19
--
※ 看板: ott 文章推薦值: 0 目前人氣: 0 累積人氣: 42
→
guest
回列表(←)
分享