以樂透為例，用Python統計馬可夫矩陣

答案是矩陣裡的每個元素都是 1/49 (下一期開出某個號碼的機率)，完畢。

寫這篇的目的並不是真的想要得到"明牌"，只是找個題目練習一下 Python 與馬可夫模型；另一方面，樂透是大家都感興趣的題目，希望這樣應該比較容易理解，並不是真的可以算出明牌，如果有執迷不悟的，本文最後會說明為什麼。

馬可夫鏈是機率問題中很常使用的模型，近年來因為大數據盛行，每個領域需要網路爬蟲收集大量資料 (分析關鍵字)，自然語言處理 (NLP) 於是成為很熱門的領域，而 Markov model 是 使用統計方法處理自然語言時常用的方法，google 也是使用 Markov model 來為每個網頁做 ranking。基本原理是，把每個網站當作節點，然後把連入與連出的連結當作是每個節點的連線，那抵達某個特定節點的機率就代表該頁面 (page) 的知名度；越容易抵達某個網頁代表被很多網頁參考，也就是該網頁知名度越高。

簡單來說，馬可夫模型主要就是從這個狀態到下一個狀態給定機率 P，而目前狀態到所有下一個狀態的機率總和一定為1。馬可夫矩陣可以用下圖表示：



上圖用大樂透來說明的話，本文要來練習的題目，目標是統計出 106 年度每個號碼跟下期號碼之間的機率，i, j 都是 49，例如: P25,49 代表本期開 25 號而下一期開 49 號的機率，而 P25,1 P25,2 ... P25,49 加總起來會是 1。

再說一次，機率理論值其實都是 1/49，這邊只是練習馬可夫模型與 Python ，因為樣本數太小，人很容易會覺得每期的號碼之間有關連，但是純粹只是機率問題。當統計的期數越多，其實會發現機率會越接近 1/49，本文只是幫助理解/說明與練習。

首先從網路上找到 106 年開獎的資料，內容如下，把內容存成 csv 檔案。Index 代表每一期開獎的號碼，Drop是落球順序，Seq是排序後的開獎順序，Special是特別號。

Index,Month,Day,Weekday,Drop1,Drop2,Drop3,Drop4,Drop5,Drop6,Seq1,Seq2,Seq3,Seq4,Seq5,Seq6,Special
1,1,3,二,18,41,47,49,22,42,18,22,41,42,47,49,45
2,,6,五,35,32,12,40,14,43,12,14,32,35,40,43,17
3,,10,二,37,14,04,46,13,39,04,13,14,37,39,46,31
4,,13,五,19,25,27,21,38,14,14,19,21,25,27,38,22
5,,17,二,40,28,12,04,34,17,04,12,17,28,34,40,07
...

我們用 read_csv() 把 csv 檔案內容載入成 dataframe，讀取從 Seq1 到 Seq6 的資料。此外，這裡假設有很多個 csv 檔案，可以使用 concat() 把兩個 dataframe 串接起來。希望兩個 dataframe 的 index 不要保留，因此設定 ignore_index 為 True。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

import numpy as np import pandas as pd if __name__ == '__main__': folder = "D:\\data" files = ['105.csv', '106.csv'] # concate and select the required csv data try: print ("Input files: " + str(files)) df_all = [] for file in files: if not file.endswith("csv"): continue path = "%s\\%s"%(folder, file) df = pd.read_csv(path) df_seq = df.loc[::1, 'Seq1':'Seq6'] if len(df_all) == 0: df_all = pd.DataFrame(columns=df.columns).loc[::1, 'Seq1':'Seq6'] df_all = pd.concat([df_all, df_seq], ignore_index=True) except Exception as e: print ("read and concate csv got exception!") raise e print (df_all)

上面這段代碼執行後會看到以下輸出：
Input files: ['106.csv']
    Seq1 Seq2 Seq3 Seq4 Seq5 Seq6
0     18   22   41   42   47   49
1     12   14   32   35   40   43
2      4   13   14   37   39   46
3     14   19   21   25   27   38
4      4   12   17   28   34   40
5      1    5   16   18   31   41
6      1    9   17   31   33   48
...

現在我們已經有了所有開獎的資料，接下來要用過去的資料統計出一個馬可夫矩陣。

第4行先用 NumPy 初始化一個 50x50 的矩陣，初始內容為 0。為了方便理解，第0行/列不會有資料。所以 stat [25, 49] 代表本期開出 25 號，下一期開出 49 號的次數。i 代表 csv 裡每一列的資料，j 是 6 個號碼其中一個。
第10行，x 代表的就是本期 (i) 的號碼，
第11行，y 就是下一期 (i+1) 的號碼。

1 2 3 4 5 6 7 8 9 10 11 12 13

# just use the numbers as index, no data in row/col 0 # x: current number, y: next number # stat[x, y] = the counts of when current index is x and the next index is y stat = np.zeros((50, 50), dtype=float) # the last data will not be used. # the last data will be used to get the furture prediction for i in range(len(df_all.values)-1): for j in range(6): x = df_all.values[i][j] for y in df_all.values[i+1]: stat[x, y] += 1.0 print(stat)

上面這段代碼會輸出如下，有需要可以再自己指定想看的陣列元素。
這邊可以看到第0行/列都不會有資料 ( 次數為0 )

[108 rows x 6 columns]
[[0. 0. 0. ... 0. 0. 0.]
 [0. 5. 1. ... 0. 2. 2.]
 [0. 2. 2. ... 2. 0. 3.]
 ...
 [0. 0. 5. ... 0. 1. 0.]
 [0. 2. 1. ... 2. 0. 1.]
 [0. 1. 0. ... 1. 1. 0.]]

接下來，若要換算成機率，如下：

1 2 3 4 5 6 7 8 9 10 11

summary = np.zeros(50, dtype=float) percent = np.zeros((50, 50), dtype=float) for i in range(50): summary[i] = np.sum(stat[i]) for j in range(50): if stat[i, j] != 0: percent[i ,j] = stat[i, j] / summary[i] print (percent) for i in range(50): print (str(i) + " max: " + str(np.max(percent[i])))

會看到以下的輸出，當樣本數越多，會看到機率越接近 1/49 (0.020...)
[[0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.06944444 0.01388889 ... 0.         0.02777778 0.02777778]
 [0.         0.02083333 0.02083333 ... 0.02083333 0.         0.03125   ]
 ...
 [0.         0.         0.07575758 ... 0.         0.01515152 0.        ]
 [0.         0.03703704 0.01851852 ... 0.03703704 0.         0.01851852]
 [0.         0.01388889 0.         ... 0.01388889 0.01388889 0.        ]]
0 max: 0.0
1 max: 0.06944444444444445
2 max: 0.0625
3 max: 0.06060606060606061
4 max: 0.047619047619047616
5 max: 0.0641025641025641
6 max: 0.06666666666666667
7 max: 0.0625
8 max: 0.052083333333333336
9 max: 0.07407407407407407
...

上面代碼我們其實省略了最後一筆的資料，不加入統計，最後，我們用最後一期的資料來找下一期對應本期號碼機率最大的五個 Pi, j。
np.argsort() 會傳回排序後元素的索引值，預設是由小排到大，所以傳回最後五個元素，就是機率最大的五個 Pi, j的索引值。( 索引值就是開獎號碼 )

1 2 3 4

last = df_all.values[len(df_all.values)-1] for i in range(len(last)): idx = np.argsort(stat[last[i]]) print ("%d : %s" % (last[i], str(idx[45:50])))

輸出如下，可以對比實際上開出的號碼。

11 : [23 33 32 11 18]
12 : [11 15 39 35 10]
15 : [28 23 41 35 16]
20 : [44 20  9 12 22]
27 : [12 49 27 18 42]
47 : [19 42 14 30  2]
correct answer: 06,12,18,36,46,47

開獎並非在實驗室進行，每次的變因很多，就算真的存在一個方程式，實際上也會因為環境因素的誤差而產生不同的結果，本文純粹只是一個練習。