※ 本文為 zbali.bbs. 轉寄自 ptt.cc 更新時間: 2020-02-06 23:26:49
看板 Gossiping
作者 Glamsight (安穩殘憶)
標題 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
 

本文就考量下一次感染的人數亦可與現有已感染人數呈正比關係

及根據 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
[圖]
 
: : 1/30
: : 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/ni4FgcQ.png
[圖]
 
: 所以明天中國大概 確診 28297 人,增加 4046 人
:                  死亡   566 人,增加   76 人
: 個人希望他們不要真的這樣玩數據。

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.112.90.228 (臺灣)
※ 文章代碼(AID): #1UF2vK-E (Gossiping)
※ 文章網址: https://www.ptt.cc/bbs/Gossiping/M.1581002324.A.F8E.html
WolfTeacher: 恩恩 跟我想的一樣1F 180.176.0.145 台灣 02/06 23:19
※ 編輯: Glamsight (140.112.90.228 臺灣), 02/06/2020 23:19:35
widec: 對!我也這麼覺得!!2F 36.232.18.241 台灣 02/06 23:19
gensokyo: 這個求證有厲害到推一個3F 140.113.170.68 台灣 02/06 23:19

--
※ 看板: ott 文章推薦值: 0 目前人氣: 0 累積人氣: 42 
作者 Glamsight 的最新發文:
點此顯示更多發文記錄
分享網址: 複製 已複製
guest
x)推文 r)回覆 e)編輯 d)刪除 M)收藏 ^x)轉錄 同主題: =)首篇 [)上篇 ])下篇