【彩票】彩票預測算法(一)：離散型馬爾可夫鏈模型C#實現

前言：彩票是一個坑，千萬不要往里面跳。任何預測彩票的方法都不可能100%，都只能說比你盲目去買要多那么一些機會而已。

　　已經3個月沒寫博客了，因為業余時間一直在研究彩票，發現還是有很多樂趣，偶爾買買，娛樂一下。本文的目的是向大家分享一個經典的數學預測算法的思路以及代碼。對於這個馬爾可夫鏈模型，我本人以前也只是聽說過，研究不深，如有錯誤，還請賜教，互相學習。

1.馬爾可夫鏈預測模型介紹

　　馬爾可夫鏈是一個能夠用數學方法就能解釋自然變化的一般規律模型，它是由著名的俄國數學家馬爾科夫在1910年左右提出的。馬爾科夫過程已經是現在概率論中隨機過程理論的一個重要方面。經過了一百年左右的發展，馬爾可夫過程已經滲透到各個領域並發揮了重要的作用，如在我們熟知的經濟、通信領域，除此之外在地質災害、醫療衛生事業、生物學等自然科學領域也發揮了非常重要的作用。

　　人們在對實際問題的研究中會發現隨着時間的持續發展變化會產生很多現象。還有一些現象或過程可以表述如下：在“現在”是已知的情況下，這種變化過程的“未來”與“過去”是毫無聯系的。也就是說這種過程的未來所出現的情況不依賴於過去的發展變化，我們就把具有上述性質的過程稱之為馬爾可夫過程。馬爾可夫過程可以描述現實生活中的很多現象。例如，我們熟知的液體中的顆粒所做的布朗運動、在商業活動中所要研究的每天銷售情況、在數字通信中的語音信號、視頻信號等。馬爾可夫鏈在其他領域的應用還有很多，如在銀行的不良資產的管理、機車管理、企業管理、生態環境演變、城市用水量仿真、信息處理等科學研究和生產生活中都有廣泛應用。

2.馬爾可夫鏈的數學概念和性質

        定義1：

　 

    定義2：

       

     上面是2個最簡單的馬爾可夫鏈的數學定義，看不懂沒關系，簡單解釋一下：

     1.從狀態k到k+1與時間k無關，也就是說這個隨機過程與時間k無關，而從k到k+1狀態，有一個轉移概率，馬爾可夫鏈的核心其實也就是這個轉移概率；

     2.根據馬爾可夫鏈的思想，一步轉移概率Pij很容易得到，但是預測的時候，往往要根據最近K期的數據來進行，所以要計算K步轉移概率；

     3.任意步的轉移概率可以根據C-K方程來計算，CK方程是一種計算轉移概率的基本方法，簡單的算法就是：通過一步轉移概率矩陣P獨自相乘m次，就可以得到m步轉移概率。

     4.馬爾可夫鏈的思想，就是根據歷史的數據，統計得到轉移概率，然后根據滯時權重對每個狀態進行預測，概率最高的是最可能出現的。

     5.對於離散型馬爾可夫鏈序列變量，一般計算之前需要對變量進行“馬氏性”檢驗，統計量就是卡方分布。

     6.馬爾可夫鏈的研究還有很多其他的方面，比如狀態分類，極限概率，平穩分布等等，這些太高級，沒時間去搞很懂，這些對預測過程的精度是有一定影響的。

3.離散型馬爾可夫鏈變量預測步驟

3.1 狀態分類

　　對於離散型變量來說，首先要把目標的數據進行歸類，對模型來說，一般狀態都是有限的，比如說雙色球，可以把16個籃球號碼分為8個狀態，2個一組。當然一些經濟和實際生活數據的狀態分類，就要根據實際情況了。

3.2 計算轉移概率矩陣

　　轉移概率矩陣是可以根據歷史數據的頻率f(i,j)統計得到。f(i,j)是狀態i到狀態j轉移的次數；然后概率轉移矩陣

                    p(i,j) = f(i,j)/f(i.)  ;頻數除以當前行的和值即為概率

3.3 "馬氏性"檢驗

　　對於離散型的變量，需要利用歷史數據進行“馬氏性”檢驗。檢驗公式為：

                                 

    然后根據顯著性水平(程序中固定取0.05) ，查表求m自由度時的閥值,若 ，則滿足 馬氏性，可以進行下一步的預測，否則沒有多大的意義。

3.4 計算自相關系數和各種步長的權重

 　若滿足馬氏性，就可以對下一個狀態進行預測了，預測根據滯時k，有權重調整，權重W(k)是根據自相關系數R(k)計算得到的，公式如下：

            

           

　　k為滯時期，我程序測試里面選的5，L是總的歷史數據次數,X是歷史數據序列。

3.5 計算不同滯時期的轉移概率矩陣

      根據C-K方程提供的算法，計算k步的轉移概率矩陣 P_i^(k) ,又一次轉移概率矩陣自乘 k次得到。

3.6 預測下一個狀態

　　下一個狀態的預測概率通過相同狀態的各個預測概率加權和得到，計算用到公式：

                       

       最后一步的時候要注意，要根據最后k期的歷史數據所在狀態值和步長的權值相乘。滯時期為1的數據，是最后1期數據(最新的數據)，這個循環的時候要注意，很容易掉進坑里。

4.離散型馬爾可夫鏈模型代碼

　　本文使用C#實現了簡單的離散型馬爾可夫鏈模型，在驗證馬氏性的時候，由於需要查表求值，所以暫時固定了自由度25,顯著性水平0.05，模型核心代碼：

 

/// <summary>離散型馬爾可夫鏈預測模型</summary>
    public class DiscreteMarkov
    {
        #region 屬性
        /// <summary>樣本點狀態時間序列,按照時間升序</summary>
        public List<int> StateList { get; set; }
        /// <summary>狀態總數,對應模型的m</summary>
        public int Count { get; set; }
        /// <summary>概率轉移矩陣Pij</summary>
        public List<DenseMatrix> ProbMatrix { get; set; }
        /// <summary>各階的自相關系數</summary>
        public double[] Rk { get; set; }
        /// <summary>各階的權重/summary>
        public double[] Wk { get; set; }
        /// <summary>頻數矩陣/summary>
        public int[][] CountStatic { get; set; }
        /// <summary>目標序列是否滿足"馬氏性"/summary>
        public Boolean IsMarkov { get; set; }
        /// <summary>滯時期，K/summary>
        public int LagPeriod { get; set; }

        /// <summary>預測概率</summary>
        public double[] PredictValue { get; set; }
        #endregion

        #region 構造函數
        public DiscreteMarkov(List<int> data, int count,int K = 5)
        {
            this.StateList = data;
            this.LagPeriod = K;
            this.Count = count;
            this.CountStatic = StaticCount(data, count);
            this.ProbMatrix = new List<DenseMatrix>();
            var t0 = DenseMatrix.OfArray(StaticProbability(this.CountStatic).ConvertToArray<double>());
            ProbMatrix.Add(t0);
          
            for (int i = 1; i < K; i++) //根據CK方程，計算各步的狀態轉移矩陣
            {
                var temp = ProbMatrix[i - 1] * t0;
                ProbMatrix.Add(temp);
            }
            if (ValidateMarkov())
            {
                CorrCoefficient();
                TimeWeight();
                PredictProb();
            }
            else
            {
                Console.WriteLine("馬氏性 檢驗失敗,無法進行下一步預測");
            }
        }
        #endregion

        #region 驗證
        /// <summary>驗證是否滿足馬氏性,默認的顯著性水平是0.05，自由度25</summary>
        /// <returns></returns>
        public Boolean ValidateMarkov()
        {
            //計算列和
            int[] cp1 = new int[Count];
            int allcount = CountStatic.Select(n => n.Sum()).Sum();//總數

            for (int i = 0; i < Count; i++)
            {
                for (int j = 0; j < Count; j++) cp1[i] += CountStatic[j][i];
            }
            double[] cp = cp1.Select(n => (double)n / (double)allcount).ToArray();

            //計算伽馬平方統計量
            double gm = 0;
            for (int i = 0; i < Count; i++)
            {
                for (int j = 0; j < Count; j++)
                {
                    if (CountStatic[i][j] != 0)
                        gm += 2 * CountStatic[i][j] * Math.Abs(Math.Log(ProbMatrix[0][i,j] / cp[j], Math.E));
                }
            }
            //查表求a = 0.05時，伽馬分布的臨界值F(m-1)^2,如果實際的gm值大於差別求得的值，則滿足
            //查表要自己做表，這里只演示0.05的情況  卡方分布            
            return gm >= 37.65;
        }

        /// <summary>計算相關系數</summary>
        public void CorrCoefficient()
        {
            double mean = (double)StateList.Sum() /(double) StateList.Count;//均值

            double p = StateList.Select(n => (n - mean)*(n-mean)).Sum();

            Rk = new double[LagPeriod];

            for (int i = 0; i < LagPeriod; i++)
            {
                double s1 = 0;              
                for (int L = 0; L < StateList.Count - LagPeriod ; L++)
                {
                    s1 += (StateList[L] - mean) * (StateList[L + i] - mean);
                }
                Rk[i] = s1 / p;
            }
        }

        /// <summary>計算滯時的步長</summary>
        public void TimeWeight()
        {
            double sum = Rk.Select(n=>Math.Abs(n)).Sum();
            Wk = Rk.Select(n => Math.Abs(n) / sum).ToArray();
        }

        /// <summary>預測狀態概率</summary>
        public void PredictProb()
        {
            PredictValue = new double[Count];
            //這里很關鍵，權重和滯時的關系要顛倒，循環計算的時候要注意
            //另外，要根據最近幾期的出現數，確定概率的狀態，必須取出最后幾期的數據

            //1.先取最后K期數據
            var last = StateList.GetRange(StateList.Count - LagPeriod, LagPeriod);
            //2.注意last數據是升序,最后一位對於的滯時期 是k =1 
            for (int i = 0; i < Count; i++)
            {
                for (int j = 0; j < LagPeriod; j++)
                {
                    //滯時期j的數據狀態
                    var state = last[last.Count - 1 - j] - 1;
                    PredictValue[i] += Wk[j] * ProbMatrix[j][state ,i];
                }
            }
        }
        #endregion

        #region 靜態 輔助方法
        /// <summary>統計頻數矩陣</summary>
        /// <param name="data">升序數據</param>
        public static int[][] StaticCount(List<int> data, int statusCount)
        {
            int[][] res = new int[statusCount][];

            for (int i = 0; i < statusCount; i++) res[i] = new int[statusCount];

            for (int i = 0; i < data.Count - 1; i++) res[data[i]-1][data[i + 1]-1]++;

            return res;
        }
        /// <summary>根據頻數，計算轉移概率矩陣</summary>
        /// <param name="data">頻率矩陣</param>
        public static double[][] StaticProbability(int[][] data)
        {
            double[][] res = new double[data.Length][];
            for (int i = 0; i < data.Length; i++)
            {
                int sum = data[i].Sum();
                res[i] = data[i].Select(n => (double)n / (double)sum).ToArray();
            }
            return res;
        }
        #endregion
    }

調用方法很簡單，如下代碼：這里使用的是論文文獻中的數據，單個號碼的隨機概率為16.6%,程序預測的概率可以到25-30%的樣子，應該還有調整的空間。

 //歷史狀態數據
            List<int> data = new List<int>(){
            6,4,4,5,2,4,6,1,2,6,  5,6,4,4,6 , 5,3,6,5,2 , 5,3,3,4,4,
            4,1,1,1,1,3,5,6,5,5,  5,5,4,6,5 , 4,1,3,1,3 , 1,3,1,2,5,
            2,2,5,5,1,4,4,2,6,1,  5,4,6,3,2,  2,6,4,4,4,  4,3,1,5,3,
            1,2,6,5,3,6,3,6,4,6,  2,4,4,6,3,  3,6,2,6,1,  3,2,2,6,6,
            4,4,3,1,4,1,2,6,4,4,  1,2};//,6,4,3,6,2,5,5,5

            var result = new DiscreteMarkov(data, 6,5);

5.實際案例

     哈哈，請關注博客，年后將根據此算法，對高頻彩快3和11選5進行實證分析。因為這個過程有點復雜，不是一下子可以搞定的。

     本文的相關文字資料，公式和數據來源根據這篇文獻：“馬爾可夫鏈預測模型及一些應用”，2012.3 溫海彬

     最后，彩票風險很大，購彩需謹慎。你的熱情和推薦，是我的動力哦。

    補充一下，其中有一個擴展方法，進行數組轉換的，忘記貼上去了：    

public static T[,] ConvertToArray<T>(this T[][] data)
        {
            T[,] res = new T[data.Length, data[0].Length];
            for (int i = 0; i < data.Length; i++)
            {
                for (int j = 0; j < data[i].Length; j++) res[i, j] = data[i][j];
            }
            return res;
        }