用Python入门不明觉厉的马尔可夫链蒙特卡罗_光环大数据python培训
马尔可夫链模型python实现
马尔可夫链模型python实现马尔可夫链模型是一种统计模型，它假设未来的状态只与当前状态有关，而与过去的状态无关。
这种特性被称为“无记忆性”或“马尔可夫性质”。
在Python中，我们可以使用NumPy和Pandas等库来实现马尔可夫链模型。
以下是一个简单的马尔可夫链模型的Python实现：pythonimport numpy as npimport pandas as pd# 假设我们有一个状态集合，状态0可以转移到状态0或状态1，状态1只能转移到状态0states = ['state0', 'state1']# 定义转移概率矩阵transition_matrix = np.array([[0.9, 0.1], # 状态0转移到状态0的概率是0.9，转移到状态1的概率是0.1[1.0, 0.0] # 状态1只能转移到状态0])# 初始化当前状态current_state = 'state0'# 模拟马尔可夫链for _ in range(10):print(f"当前状态: {current_state}")# 根据转移概率矩阵确定下一个状态next_state_probabilities =transition_matrix[states.index(current_state)]next_state = np.random.choice(states, p=next_state_probabilities)current_state = next_state# 输出最终的状态print(f"最终状态: {current_state}")上述代码首先定义了状态集合和转移概率矩阵，然后初始化了当前状态。
接着，它使用一个循环来模拟马尔可夫链的演化，每次迭代都会根据当前的状态和转移概率矩阵来确定下一个状态。
最后，它输出了最终的状态。
马尔可夫链模型python实现
马尔可夫链模型python实现全文共四篇示例，供读者参考第一篇示例：马尔可夫链是一种随机过程，它基于马尔可夫性质，即未来的状态只取决于当前的状态，而不受过去的影响。
马尔可夫链模型广泛应用于自然语言处理、机器学习、统计建模等领域，可以用来模拟具有随机性的现象。
在本文中，我们将介绍如何使用Python实现马尔可夫链模型。
我们需要了解马尔可夫链的基本概念。
马尔可夫链由状态空间、初始状态和状态转移概率矩阵组成。
状态空间是所有可能状态的集合，初始状态指定了链条起始状态，状态转移概率矩阵描述了从一个状态到另一个状态的转移概率。
接下来，我们将通过一个简单的例子来说明如何使用Python实现马尔可夫链模型。
假设我们有一个天气预测的问题，天气状态包括“晴天”和“雨天”，我们希望根据过去的天气情况预测未来的天气。
我们需要定义状态空间和状态转移概率矩阵。
状态空间定义如下：接着，我们可以定义状态转移概率矩阵，假设转移概率如下：以上代码中的transition_matrix表示在晴天时，下一天为晴天的概率为0.8，为雨天的概率为0.2；在雨天时，下一天为晴天的概率为0.4，为雨天的概率为0.6。
接着，我们可以编写Python代码来实现马尔可夫链模型。
我们需要定义一个函数来根据当前状态和转移概率矩阵来确定下一个状态：```pythonimport randomdef next_state(current_state, transition_matrix):next_states = transition_matrix[current_state]probabilities = list(next_states.values())next_state = random.choices(list(next_states.keys()), weights=probabilities)[0]return next_state```以上代码定义了一个next_state函数，接受当前状态和转移概率矩阵作为参数，返回根据转移概率确定的下一个状态。
马尔可夫链蒙特卡洛方法简介(七)
马尔可夫链蒙特卡洛方法简介蒙特卡洛方法是一种通过随机抽样来解决问题的数值计算方法。
而在蒙特卡洛方法中，马尔可夫链蒙特卡洛方法（Markov Chain Monte Carlo, MCMC）是一种重要的技术，它可以用于求解很多实际问题，比如概率分布的估计、贝叶斯统计推断等。
本文将对马尔可夫链蒙特卡洛方法进行简要介绍。
1. 马尔可夫链马尔可夫链是指一个具有马尔可夫性质的随机过程。
所谓马尔可夫性质是指一个系统在给定当前状态下，未来的状态只与当前状态有关，而与过去状态无关。
换句话说，马尔可夫链的未来状态只取决于当前状态，而与过去状态无关。
这种性质使得马尔可夫链在模拟复杂系统时非常有用。
2. 马尔可夫链蒙特卡洛方法在蒙特卡洛方法中，马尔可夫链蒙特卡洛方法是通过构造一个马尔可夫链，使得该链的平稳分布恰好是我们要求的概率分布。
通过对该马尔可夫链进行随机抽样，最终可以得到与平稳分布一致的样本，从而对概率分布进行估计。
3. Metropolis-Hastings算法Metropolis-Hastings算法是一种常用的马尔可夫链蒙特卡洛方法。
其基本思想是通过一系列状态转移来构造一个满足平稳分布的马尔可夫链。
具体而言，算法首先随机初始化一个状态，然后通过一定的转移规则来进行状态转移。
在每次状态转移后，我们都根据一定的准则来接受或者拒绝转移，以保证最终的样本满足平稳分布。
4. Gibbs采样Gibbs采样是一种特殊的Metropolis-Hastings算法。
它适用于高维参数的分布估计问题。
在Gibbs采样中，我们将多维参数分解为多个条件分布，然后通过依次对每个条件分布进行抽样来得到最终的样本。
Gibbs采样在贝叶斯统计推断等领域有着广泛的应用。
5. 贝叶斯统计推断马尔可夫链蒙特卡洛方法在贝叶斯统计推断中有着重要的应用。
在贝叶斯统计中，我们往往需要对参数的后验分布进行估计。
而马尔可夫链蒙特卡洛方法可以通过对后验分布进行抽样来进行估计，从而得到参数的后验分布的近似值。
小白都能看懂的蒙特卡洛方法以及python实现
小白都能看懂的蒙特卡洛方法以及python实现蒙特卡洛方法是一种基于随机抽样的数学方法，它通过模拟随机过程来求解复杂问题。
这种方法在许多领域都有应用，例如计算机科学、物理学、金融学等。
今天，我们就来介绍一下小白都能看懂的蒙特卡洛方法以及python实现。
一、什么是蒙特卡洛方法？蒙特卡洛方法是一种基于概率的数值分析方法，它通过模拟随机过程来求解复杂问题。
这种方法的基本思想是通过随机抽样来估计一个未知量的数值。
在蒙特卡洛方法中，我们通常会建立一个概率模型，模拟随机过程，并通过对模型进行大量的抽样，来估计未知量的数值。
二、为什么要用蒙特卡洛方法？蒙特卡洛方法具有许多优点，例如计算速度快、适用范围广、易于实现等。
在许多实际问题中，我们无法直接求解数学模型，而蒙特卡洛方法可以通过模拟随机过程来求解复杂问题，从而得到近似解。
此外，蒙特卡洛方法还可以用于解决一些难以用传统数学方法解决的问题。
三、Python实现蒙特卡洛方法下面是一个简单的Python代码示例，演示了如何使用蒙特卡洛方法估算圆周率π的值：```pythonimportrandomimportmathdefestimate_pi(n):#创建一个正方形区域，并随机生成点在区域内points=[(random.uniform(0,1),random.uniform(0,1))for_inra nge(n)]#将点落在正方形区域内的圆心角缩小到π/n弧度内foriinrange(n):x,y=points[i]dx,dy=x*2,y*2points[i]=(x+dx*math.sin(math.pi/n*(i+1)),y+dy*math.cos(m ath.pi/n*(i+1)))#统计落在圆内的点数inside_points=len([pforpinpointsifmath.sqrt(math.pow(p[0] -0,2)+math.pow(p[1]-0,2))<=1])#估算π的值pi_estimate=4*inside_points/nreturnpi_estimate```这段代码中，我们首先创建了一个正方形区域，并随机生成了一些点在区域内。
从零开始学Python_光环大数据python培训
从零开始学Python_光环大数据python培训一、频数统计我们以被调查用户的收入数据为例，来谈谈频数统计函数value_counts。
频数统计，顾名思义就是统计某个离散变量各水平的频次。
这里统计的是性别男女的人数，是一个绝对值，如果想进一步查看男女的百分比例，可以通过下面的方式实现：而在R语言中，table函数就是起到频数统计的作用，另外还提供了更加灵活的prop.table函数，可以直接求出比例。
如上是单变量的频数统计，如果需要统计两个离散变量的交叉统计表，该如何实现？不急，pandas模块提供了crosstab函数，我们来看看其用法：R语言的话，任然使用table函数即可。
二、缺失值处理在数据分析或建模过程中，我们希望数据集是干净的，没有缺失、异常之类，但面临的实际情况确实数据集很脏，例如对于缺失值我们该如何解决？一般情况，缺失值可以通过删除或替补的方式来处理。
首先是要监控每个变量是否存在缺失，缺失的比例如何？这里我们借助于pandas模块中的isnull函数、dropna函数和fillna函数。
首先，我们手工编造一个含缺失值的数据框：其次，使用isnull函数检查数据集的缺失情况：最后，对缺失数据进行处理：删除法dropna函数，有两种删除模式，一种是对含有缺失的行(任意一列)进行删除，另一种是删除那些全是缺失(所有列)的行，具体如下：由于df数据集不存在行全为缺失的观测，故没有实现删除。
替补法fillna函数提供前向替补、后向替补和函数替补的几种方法，具体可参见下面的代码示例：再来看看R语言是如何重现上面的操作的：不幸的是，R中没有删除每行元素都是缺失的观测，我们自定义个函数也可以实现：关于缺失值的替补，在R语言中可以使用Hmisc包中的impute函数，具体操作如下：三、数据映射大家都知道，Python和R在做循环时，效率还是很低的，如何避开循环达到相同的效果呢？这就是接下来我们要研究的映射函数apply。
从零开始学Python_光环大数据分析培训
从零开始学Python_光环大数据分析培训使用numpy构建矩阵数组的创建可以使用numpy模块中的array函数实现，一维数组只需要给array函数传入一个列表或元组，二维数组则是传入嵌套的列表或元组。
具体举例可知：arr1和arr2为一维数组，arr3为二维数组，返回一个数组的行数和列数可使用shape方法，即元素的获取使用索引的方式，查询一维数组和二维数组的元素。
一维数组的索引与列表、元组的索引完全一致，这里就不在赘述；二维数组的索引就稍微有点复杂，我们可以通过例子来说明：print函数中的‘/n’，目的用来换行，使打印出来的结果不显得那么拥挤。
咦？报告，你最后一个返回的结果错了，你不是要返回由第一行、第三行、第三列和第四列组成的2×2矩阵吗？为什么是一个1×2的一维数组？如果像上面红框中使用索引的话，将获取【0,2】和【2,3】对应的两个值。
那该如何返回想要的2×2的矩阵呢？我们可以这样写：数学函数# 取绝对值np.absnp.fabs# 算术平方根np.sqrt# 平方np.square# 指数np.exp# 对数np.log2np.log10np.log(x,base)# 符号函数（大于0的数返回1、小于0的数返回-1、0返回0值）np.sign # 向上取整np.cell# 向下取整np.floor# 返回最近的整数np.rint# 判断是否缺失np.isnan# 判断是否有限np.isfinite# 判断是否无限np.isinf# 幂运算np.power # 余数np.mod统计函数# 最大值np.max# 浮点型的最大值np.fmax# 最小值np.mim# 浮点型的最小值np.fmin# 求和np.sum# 均值np.mean# 标准差np.std# 方差np.var# 中位数np.median映射函数apply_along_axisapply_along_axis函数与R语言中的apply函数用法一致，可以针对某个轴的方向进行函数操作，同样，而且在pandas模块中的DataFrmae对象中，可以使用apply函数达到相同的效果。
Python科学计算 - Numpy快速入门 光环大数据
Python科学计算- Numpy快速入门光环大数据光环大数据Python培训了解到，Numpy是Python的一个科学计算的库，提供了矩阵运算的功能，其一般与Scipy、matplotlib一起使用。
它可用来存储和处理大型矩阵，比Python自身的嵌套列表（nestedliststructure)结构要高效的多（该结构也可以用来表示矩阵（matrix））。
NumPy（NumericPython）提供了许多高级的数值编程工具，如：矩阵数据类型、矢量处理，以及精密的运算库。
专为进行严格的数字处理而产生。
多为很多大型金融公司使用，以及核心的科学计算组织如：LawrenceLivermore，NASA用其处理一些本来使用C++，Fortran或Matlab等所做的任务。
多维数组多维数组的类型是：numpy.ndarray使用numpy.array方法以list或tuple变量为参数产生一维数组：>>>print(np.array([1,2,3,4]))[1234]>>>print(np.array((1.2,2,3,4)) )[1.22.3.4.]>>>printtype(np.array((1.2,2,3,4)))<type'numpy.ndarray'> 以list或tuple变量为元素产生二维数组：>>>print(np.array([[1,2],[3,4]]))[[12][34]]指定数据类型例如numpy.int32,numpy.int16,andnumpy.float64等：>>>printnp.array((1.2,2,3,4),dtype=np.int32)[1234]使用numpy.arange方法>>>print(np.arange(15))[01234567891011121314]>>>printtype(np.aran ge(15))<type'numpy.ndarray'>>>>printnp.arange(15).reshape(3,5)[[01234 ][56789][1011121314]]>>>printtype(np.arange(15).reshape(3,5))<type'nu mpy.ndarray'>使用numpy.linspace方法例如，在从1到3中产生9个数：>>>print(np.linspace(1,3,10))[1.1.222222221.444444441.666666671.8 88888892.111111112.333333332.555555562.777777783.]构造特定的矩阵使用numpy.zeros，numpy.ones，numpy.eye可以构造特定的矩阵>>>print(np.zeros((3,4)))[[0.0.0.0.][0.0.0.0.][0.0.0.0.]]>>>print (np.ones((4,3)))[[1.1.1.][1.1.1.][1.1.1.][1.1.1.]]>>>print(np.eye(4)) [[1.0.0.0.][0.1.0.0.][0.0.1.0.][0.0.0.1.]]创建一个三维数组：>>>print(np.ones((3,3,3)))[[[1.1.1.][1.1.1.][1.1.1.]][[1.1.1.][1.1.1.][1.1.1.]][[1.1.1.][1.1.1.][1.1.1.]]]获取数组的属性>>>a=np.zeros((2,3,2))>>>print(a.ndim)#数组的维数3>>>print(a.shape)#数组每一维的大小(2,3,2)>>>print(a.size)#数组的元素数12>>>print(a.dtype)#元素类型float64>>>print(a.itemsize)#每个元素所占的字节数8数组索引，切片，赋值>>>a=np.array([[2,3,4],[5,6,7]])>>>print(a)[[234][567]]>>>print(a [1,2])#index从0开始7>>>printa[1,:][567]>>>print(a[1,1:2])[6]>>>a[1,:]=[8,9,10]#直接赋值>>>print(a)[[234][8910]]使用for操作元素>>>forxinnp.linspace(1,3,3):...print(x)...1.02.03.0基本的数组运算先构造数组a、b：>>>a=np.ones((2,2))>>>b=np.eye(2)>>>print(a)[[1.1.][1.1.]]>>>prin t(b)[[1.0.][0.1.]]数组的加减乘除>>>print(a>2)[[FalseFalse][FalseFalse]]>>>print(a+b)[[2.1.][1.2.] ]>>>print(a-b)[[0.1.][1.0.]]>>>print(b*2)[[2.0.][0.2.]]>>>print((a*2) *(b*2))[[4.0.][0.4.]]>>>print(b/(a*2))[[0.50.][0.0.5]]>>>print((b*2)* *4)[[16.0][016.]]使用数组对象自带的方法>>>a.sum()#a的元素个数4.0>>>a.sum(axis=0)#计算每一列（二维数组中类似于矩阵的列）的和array([2.,2.])>>>a.min()1.0>>>a.max()1.0使用numpy 下的方法>>>np.sin(a)array([[0.84147098,0.84147098],[0.84147098,0.84147098]] )>>>np.max(a)1.0>>>np.floor(a)array([[1.,1.],[1.,1.]])>>>np.exp(a)arr ay([[2.71828183,2.71828183],[2.71828183,2.71828183]])>>>np.dot(a,a)##矩阵乘法array([[2.,2.],[2.,2.]])合并数组使用numpy下的vstack和hstack函数：>>>a=np.ones((2,2))>>>b=np.eye(2)>>>print(np.vstack((a,b)))#顾名思义v--vertical垂直[[1.1.][1.1.][1.0.][0.1.]]>>>print(np.hstack((a,b)))#顾名思义h--horizonal水平[[1.1.1.0.][1.1.0.1.]]看一下这两个函数有没有涉及到浅拷贝这种问题：>>>c=np.hstack((a,b))>>>printc[[1.1.1.0.][1.1.0.1.]]>>>a[1,1]=5>> >b[1,1]=5>>>printc[[1.1.1.0.][1.1.0.1.]]可以看到，a、b中元素的改变并未影响c。
蒙特卡洛模拟python代码
蒙特卡洛模拟python代码蒙特卡罗模拟是个有效的研究工具，可以模拟复杂的系统进行研究，它通常用于收集系统的信息，帮助做出合理的决策。
在用python实现蒙特卡洛模拟时，首先要创建模拟过程所需的环境，这需要从现实环境中发现模型，并将其转化成可用的数学模型，如给定的代码：from random import random # 随机函数import matplotlib.pyplot as plt # 画图函数# 模拟的蒙特卡罗过程simulationSize = 10000 # 模拟步数inCircle = 0 # 在圆内的模拟步for i in range(simulationSize):x = random() # 产生0-1之间的随机数y = random()if x*x + y*y <= 1:inCircle += 1piValue = 4 * (inCircle / simulationSize) # 获取pi的模拟值# 绘出结果print(piValue)plt.title('Pi value')plt.plot(piValue, 'bo')plt.show()从上述代码可以看出，首先，将随机变量要求转换为均匀分布，以产生0到1之间的随机数，然后模拟仿真过程并判断点是否位于圆形内，最后，获得pi的模拟值，并使用matplotlib将模拟结果绘出。
蒙特卡洛模拟是一种重要而有效的研究工具，它在有限模拟步数，大量随机数和有效算法辅助下，可以推测复杂模型和复杂系统的运行，并且可以得出较为可靠的模拟结果。
用python实现蒙特卡罗模拟，从模型发现到技术实施，可以最大限度地利用python语言的特性，实现简洁、优雅的编程，以提高模拟的效率及正确率。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
用Python入门不明觉厉的马尔可夫链蒙特卡罗_光环大数据python培训在过去几个月里，我在数据科学的世界里反复遇到一个词：马尔可夫链蒙特卡洛（Markov Chain Monte Carlo , MCMC）。
在我的研究室、podcast和文章里，每每遇到这个词我都会“不明觉厉”地点点头，觉得这个算法听起来很酷，但每次听人提起也只是有个模模糊糊的概念。
我屡次尝试学习MCMC和贝叶斯推论，而一拿起书，又很快就放弃了。
无奈之下，我选择了学习任何新东西最佳的方法：应用到一个实际问题中。
通过使用一些我曾试图分析的睡眠数据和一本实操类的、基于应用教学的书（《写给开发者的贝叶斯方法》，我最终通过一个实际项目搞明白了MCMC。
《写给开发者的贝叶斯方法》和学习其他东西一样，当我把这些技术性的概念应用于一个实际问题中而不是单纯地通过看书去了解这些抽象概念，我更容易理解这些知识，并且更享受学习的过程。
这篇文章介绍了马尔可夫链蒙特卡洛在Python中入门级的应用操作，这个实际应用最终也使我学会使用这个强大的建模分析工具。
此项目全部的代码和数据：https:///WillKoehrsen/ai-projects/blob/master/bayesian/ bayesian_inference.ipynb这篇文章侧重于应用和结果，因此很多知识点只会粗浅的介绍，但对于那些想了解更多知识的读者，在文章也尝试提供了一些学习链接。
案例简介我的智能手环在我入睡和起床时会根据心率和运动进行记录。
它不是100%准确的，但现实世界中的数据永远不可能是完美的，不过我们依然可以运用正确的模型从这些噪声数据提取出有价值的信息。
典型睡眠数据这个项目的目的在于运用睡眠数据建立一个能够确立睡眠相对于时间的后验分布模型。
由于时间是个连续变量，我们无法知道后验分布的具体表达式，因此我们转向能够近似后验分布的算法，比如马尔可夫链蒙特卡洛（MCMC）。
选择一个概率分布在我们开始MCMC之前，我们需要为睡眠的后验分布模型选择一个合适的函数。
一种简单的做法是观察数据所呈现的图像。
下图呈现了当我入睡时时间函数的数据分布。
睡眠数据每个数据点用一个点来表示，点的密度展现了在固定时刻的观测个数。
我的智能手表只记录我入睡的那个时刻，因此为了扩展数据，我在每分钟的两端添加了数据点。
如果我的手表记录我在晚上10：05入睡，那么所有在此之前的时间点被记为0（醒着），所有在此之后的时间点记为1（睡着）。
这样一来，原本大约60夜的观察量被扩展成11340个数据点。
可以看到我趋向于在10：00后几分钟入睡，但我们希望建立一个把从醒到入睡的转变用概率进行表达的模型。
我们可以用一个简单的阶梯函数作为模型，在一个精确时间点从醒着（0）变到入睡（1），但这不能代表数据中的不确定性。
我不会每天在同一时间入睡，因此我们需要一个能够模拟出这个个渐变过程的函数来展现变化当中的差异性。
在现有数据下最好的选择是logistic函数，在0到1之前平滑地移动。
下面这个公式是睡眠状态相对时间的概率分布，也是一个logistic公式。
在这里，β (beta) 和α (alpha) 是模型的参数，我们只能通过MCMC去模拟它们的值。
下面展示了一个参数变化的logistic函数。
一个logistic函数能够很好的拟合数据，因为在logistic函数中入睡的概率在逐渐改变，捕捉了我睡眠模式的变化性。
我们希望能够带入一个具体的时间t到函数中，从而得到一个在0到1之间的睡眠状态的概率分布。
我们并不会直接得到我是否在10：00睡着了的准确答案，而是一个概率。
创建这个模型，我们通过数据和马尔可夫链蒙特卡洛去寻找最优的alpha和beta系数估计。
马尔可夫链蒙特卡洛马尔可夫链蒙特卡罗是一组从概率分布中抽样，从而建立最近似原分布的函数的方法。
因为我们不能直接去计算logistic分布，所以我们为系数(alpha 和beta)生成成千上万的数值-被称为样本-去建立分布的一个模拟。
MCMC背后的基本思想就是当我们生成越多的样本，我们的模拟就更近似于真实的分布。
马尔可夫链蒙特卡洛由两部分组成。
蒙特卡洛代表运用重复随机的样本来获取一个准确答案的一种模拟方法。
蒙特卡洛可以被看做大量重复性的实验，每次更改变量的值并观察结果。
通过选择随机的数值，我们可以在系数的范围空间，也就是变量可取值的范围，更大比例地探索。
下图展示了在我们的问题中，一个使用高斯分布作为先验的系数空间。
能够清楚地看到我们不能在这些图中一一找出单个的点，但通过在更高概率的区域（红色）进行随机抽样，我们就能够建立最近似的模型。
马尔可夫链（Markov Chain）马尔可夫链是一个“下个状态值只取决于当前状态”的过程。
（在这里，一个状态指代当前时间系数的数值分配）。
一个马尔可夫链是“健忘”的，因为如何到达当前状态并不要紧，只有当前的状态值是关键。
如果这有些难以理解的话，让我们来设想一个每天都会经历的情景–天气。
如果我们希望预测明天的天气，那么仅仅使用今天的天气状况我们就能够得到一个较为合理的预测。
如果今天下雪，我们可以观测有关下雪后第二天天气的历史数据去预测明天各种天气状况的可能性。
马尔可夫链的定义就是我们不需要知道一个过程中的全部历史状态去预测下一节点的状态，这种近似在许多现实问题中都很有用。
把马尔可夫链（Markov Chain）和蒙特卡洛（Monte Carlo），两者放到一起，就有了MCMC。
MCMC是一种基于当前值，重复为概率分布系数抽取随机数值的方法。
每个样本都是随机的，但是数值的选择也受当前值和系数先验分布的影响。
MCMC可以被看做是一个最终趋于真实分布的随机游走。
为了能够抽取alpha 和 beta的随机值，我们需要为每个系数假设一个先验分布。
因为我们没有对于这两个系数的任何假设，我们可以使用正太分布作为先验。
正太分布，也称高斯分布，是由均值（展示数据分布），和方差（展示离散性）来定义的。
下图展示了多个不同均值和离散型的正态分布。
具体的MCMC算法被称作Metropolis Hastings。
为了连接我们的观察数据到模型中，每次一组随机值被抽取，算法将把它们与数据进行比较。
一旦它们与数据不吻合（在这里我简化了一部分内容），这些值就会被舍弃，模型将停留在当前的状态值。
如果这些随机值与数据吻合，那么这些值就被接纳为各个系数新的值，成为当前的状态值。
这个过程会有一个提前设置好的迭代次数，次数越多，模型的精确度也就越高。
把上面介绍的整合到一起，就能得到在我们的问题中所需进行的最基本的MCMC步骤：为logistic函数的系数alpha 和beta选择初始值。
基于当前状态值随机分配给alpha 和beta新的值。
检查新的随机值是否与观察数据吻合。
如果不是，舍弃掉这个值，并回到上一状态值。
如果吻合，接受这个新的值作为当前状态值。
重复步骤2和3（重复次数提前设置好）。
这个算法会给出所有它所生成的alpha 和beta值。
我们可以用这些值的平均数作为alpha 和beta在logistic函数中可能性最大的终值。
MCMC不会返回“真实”的数值，而是函数分布的近似值。
睡眠状态概率分布的最终模型将会是以alph和beta均值作为系数的logistic函数。
Python实施我再三思考模拟上面提到的细节，最终我开始用Python将它们变成现实。
观察一手的结果会比阅读别人的经验贴有帮助得多。
想要在Python中实施MCMC，我们需要用到PyMC3贝叶斯库，它省略了很多细节，方便我们创建模型，避免迷失在理论之中。
通过下面的这些代码可以创建完整的模型，其中包含了参数alpha 、beta、概率p以及观测值observed，step变量是指特定的算法，sleep_trace包含了模型创建的所有参数值。
with pm.Model() as sleep_model: # Create the alpha and beta parameters # Assume a normal distribution alpha = pm.Normal(‘alpha’, mu=0.0, tau=0.05, testval=0.0) beta = pm.Normal(‘beta’, mu=0.0, tau=0.05, testval=0.0) # The sleep probability is modeled as a logistic function p = pm.Deterministic(‘p’,1./ (1.+ tt.exp(beta * time + alpha))) # Create the bernoulli parameter which uses observed data to inform the algorithm observed = pm.Bernoulli(‘obs’, p, observed=sleep_obs) # Using Metropolis Hastings Sampling step = pm.Metropolis() # Draw the specified number of samples sleep_trace = pm.sample(N_SAMPLES, step=step);为了更直观地展现代码运行的效果，我们可以看一下模型运行时alpha和beta生成的值。
这些图叫做轨迹图，可以看到每个状态都与其历史状态相关，即马尔可夫链；同时每个值剧烈波动，即蒙特卡洛抽样。
使用MCMC时，常常需要放弃轨迹图中90%的值。
这个算法并不能立即展现真实的分布情况，最初生成的值往往是不准确的。
随着算法的运行，后产生的参数值才是我们真正需要用来建模的值。
我使用了一万个样本，放弃了前50%的值，但真正在行业中应用时，样本量可达成千上万个、甚至上百万个。
通过足够多的迭代，MCMC逐渐趋近于真实的值，但是估算收敛性并不容易。
这篇文章中并不会涉及到具体的估算方法（方法之一就是计算轨迹的自我相关性），但是这是得到最准确结果的必要条件。
PyMC3的函数能够评估模型的质量，包括对轨迹、自相关图的评估。
pm.traceplot(sleep_trace,[‘alpha’,’beta’])pm.autocorrplot(sleep_trace,[‘alpha’,’beta’])轨迹图（左）和自相关性图（右）睡眠模型建模、模型运行完成后，该最终结果上场了。
我们将使用最终的5000个alpha和beta值作为参数的可能值，最终创建了一个单曲线来展现事后睡眠概率：基于5000个样本的睡眠概率分布这个模型能够很好的代表样本数据，并且展现了我睡眠模式的内在变异性。