已经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步的转移概率矩阵 Pi(k) ,又一次转移概率矩阵自乘 k次得到。
3.6 预测下一个状态
下一个状态的预测概率通过相同状态的各个预测概率加权和得到，计算用到公式：
最后一步的时候要注意，要根据最后k期的历史数据所在状态值和步长的权值相乘。滞时期为1的数据，是最后1期数据(最新的数据)，这个循环的时候要注意，很容易掉进坑里。
4.离散型马尔可夫链模型代码
本文使用C#实现了简单的离散型马尔可夫链模型，在验证马氏性的时候，由于需要查表求值，所以暂时固定了*度25,显著性水平0.05，模型核心代码：
1 /// <summary>离散型马尔可夫链预测模型</summary> 2 public class DiscreteMarkov 3 { 4 #region 属性 5 /// <summary>样本点状态时间序列,按照时间升序</summary> 6 public List<int> StateList { get; set; } 7 /// <summary>状态总数,对应模型的m</summary> 8 public int Count { get; set; } 9 /// <summary>概率转移矩阵Pij</summary> 10 public List<DenseMatrix> ProbMatrix { get; set; } 11 /// <summary>各阶的自相关系数</summary> 12 public double[] Rk { get; set; } 13 /// <summary>各阶的权重/summary> 14 public double[] Wk { get; set; } 15 /// <summary>频数矩阵/summary> 16 public int[][] CountStatic { get; set; } 17 /// <summary>目标序列是否满足"马氏性"/summary> 18 public Boolean IsMarkov { get; set; } 19 /// <summary>滞时期，K/summary> 20 public int LagPeriod { get; set; } 21 22 /// <summary>预测概率</summary> 23 public double[] PredictValue { get; set; } 24 #endregion 25 26 #region 构造函数 27 public DiscreteMarkov(List<int> data, int count,int K = 5) 28 { 29 this.StateList = data; 30 this.LagPeriod = K; 31 this.Count = count; 32 this.CountStatic = StaticCount(data, count); 33 this.ProbMatrix = new List<DenseMatrix>(); 34 var t0 = DenseMatrix.OfArray(StaticProbability(this.CountStatic).ConvertToArray<double>()); 35 ProbMatrix.Add(t0); 36 37 for (int i = 1; i < K; i++) //根据CK方程，计算各步的状态转移矩阵 38 { 39 var temp = ProbMatrix[i - 1] * t0; 40 ProbMatrix.Add(temp); 41 } 42 if (ValidateMarkov()) 43 { 44 CorrCoefficient(); 45 TimeWeight(); 46 PredictProb(); 47 } 48 else 49 { 50 Console.WriteLine("马氏性 检验失败,无法进行下一步预测"); 51 } 52 } 53 #endregion 54 55 #region 验证 56 /// <summary>验证是否满足马氏性,默认的显著性水平是0.05，*度25</summary> 57 /// <returns></returns> 58 public Boolean ValidateMarkov() 59 { 60 //计算列和 61 int[] cp1 = new int[Count]; 62 int allcount = CountStatic.Select(n => n.Sum()).Sum();//总数 63 64 for (int i = 0; i < Count; i++) 65 { 66 for (int j = 0; j < Count; j++) cp1[i] += CountStatic[j][i]; 67 } 68 double[] cp = cp1.Select(n => (double)n / (double)allcount).ToArray(); 69 70 //计算伽马平方统计量 71 double gm = 0; 72 for (int i = 0; i < Count; i++) 73 { 74 for (int j = 0; j < Count; j++) 75 { 76 if (CountStatic[i][j] != 0) 77 gm += 2 * CountStatic[i][j] * Math.Abs(Math.Log(ProbMatrix[0][i,j] / cp[j], Math.E)); 78 } 79 } 80 //查表求a = 0.05时，伽马分布的临界值F(m-1)^2,如果实际的gm值大于差别求得的值，则满足 81 //查表要自己做表，这里只演示0.05的情况 卡方分布 82 return gm >= 37.65; 83 } 84 85 /// <summary>计算相关系数</summary> 86 public void CorrCoefficient() 87 { 88 double mean = (double)StateList.Sum() /(double) StateList.Count;//均值 89 90 double p = StateList.Select(n => (n - mean)*(n-mean)).Sum(); 91 92 Rk = new double[LagPeriod]; 93 94 for (int i = 0; i < LagPeriod; i++) 95 { 96 double s1 = 0; 97 for (int L = 0; L < StateList.Count - LagPeriod ; L++) 98 { 99 s1 += (StateList[L] - mean) * (StateList[L + i] - mean); 100 } 101 Rk[i] = s1 / p; 102 } 103 } 104 105 /// <summary>计算滞时的步长</summary> 106 public void TimeWeight() 107 { 108 double sum = Rk.Select(n=>Math.Abs(n)).Sum(); 109 Wk = Rk.Select(n => Math.Abs(n) / sum).ToArray(); 110 } 111 112 /// <summary>预测状态概率</summary> 113 public void PredictProb() 114 { 115 PredictValue = new double[Count]; 116 //这里很关键，权重和滞时的关系要颠倒，循环计算的时候要注意 117 //另外，要根据最近几期的出现数，确定概率的状态，必须取出最后几期的数据 118 119 //1.先取最后K期数据 120 var last = StateList.GetRange(StateList.Count - LagPeriod, LagPeriod); 121 //2.注意last数据是升序,最后一位对于的滞时期 是k =1 122 for (int i = 0; i < Count; i++) 123 { 124 for (int j = 0; j < LagPeriod; j++) 125 { 126 //滞时期j的数据状态 127 var state = last[last.Count - 1 - j] - 1; 128 PredictValue[i] += Wk[j] * ProbMatrix[j][state ,i]; 129 } 130 } 131 } 132 #endregion 133 134 #region 静态 辅助方法 135 /// <summary>统计频数矩阵</summary> 136 /// <param name="data">升序数据</param> 137 public static int[][] StaticCount(List<int> data, int statusCount) 138 { 139 int[][] res = new int[statusCount][]; 140 141 for (int i = 0; i < statusCount; i++) res[i] = new int[statusCount]; 142 143 for (int i = 0; i < data.Count - 1; i++) res[data[i]-1][data[i + 1]-1]++; 144 145 return res; 146 } 147 /// <summary>根据频数，计算转移概率矩阵</summary> 148 /// <param name="data">频率矩阵</param> 149 public static double[][] StaticProbability(int[][] data) 150 { 151 double[][] res = new double[data.Length][]; 152 for (int i = 0; i < data.Length; i++) 153 { 154 int sum = data[i].Sum(); 155 res[i] = data[i].Select(n => (double)n / (double)sum).ToArray(); 156 } 157 return res; 158 } 159 #endregion 160 }
调用方法很简单，如下代码：这里使用的是论文文献中的数据，单个号码的随机概率为16.6%,程序预测的概率可以到25-30%的样子，应该还有调整的空间。
1 //历史状态数据 2 List<int> data = new List<int>(){ 3 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 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, 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, 6 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, 7 4,4,3,1,4,1,2,6,4,4, 1,2};//,6,4,3,6,2,5,5,5 8 9 var result = new DiscreteMarkov(data, 6,5);
5.实际案例
哈哈，请关注博客，年后将根据此算法，对高频彩快3和11选5进行实证分析。因为这个过程有点复杂，不是一下子可以搞定的。
本文的相关文字资料，公式和数据来源根据这篇文献：“马尔可夫链预测模型及一些应用”，2012.3 温海彬
文中的矩阵相乘用到了Math.NET组件。这里下载：http://www.kuaipan.cn/file/id_4279479643945457.htm?source=1
有玩高频彩(快3，11选5)的朋友，可以到群（群号：280182639），一起交流心得。
最后，彩票风险很大，购彩需谨慎。你的热情和推荐，是我的动力哦。
补充一下，其中有一个扩展方法，进行数组转换的，忘记贴上去了：
1 public static T[,] ConvertToArray<T>(this T[][] data) 2 { 3 T[,] res = new T[data.Length, data[0].Length]; 4 for (int i = 0; i < data.Length; i++) 5 { 6 for (int j = 0; j < data[i].Length; j++) res[i, j] = data[i][j]; 7 } 8 return res; 9 }