问题类型1：参数估计

真实值是否等于X？

给出数据，对于参数，可能的值的概率分布是多少？

例子1：抛硬币问题

硬币扔了n次，正面朝上是h次。

参数问题

想知道 p 的可能性。给定 n 扔的次数和 h 正面朝上次数，p 的值很可能接近 0.5，比如说在 [0.48，0.52]？

说明

参数的先验信念：p∼Uniform(0,1)
似然函数：data∼Bernoulli(p)

import pymc3 as pm
import numpy.random as npr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from collections import Counter
import seaborn as snssns.set_style('white')
sns.set_context('poster')%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import warnings
warnings.filterwarnings('ignore')
from random import shuffle
total = 30
n_heads = 11
n_tails = total - n_heads
tosses = [1] * n_heads + [0] * n_tails
shuffle(tosses)

可视化数据：

def plot_coins():fig = plt.figure()ax = fig.add_subplot(1,1,1)ax.bar(list(Counter(tosses).keys()), list(Counter(tosses).values()))ax.set_xticks([0, 1])ax.set_xticklabels(['tails', 'heads'])ax.set_ylim(0, 20)ax.set_yticks(np.arange(0, 21, 5))        return figfig = plot_coins()
plt.show()

建立模型：

with pm.Model() as coin_model:# Specify prior using Uniform object.p_prior = pm.Uniform('p', 0, 1)  # Specify likelihood using Bernoulli object.like = pm.Bernoulli('likelihood', p=p_prior, observed=tosses)     # "observed=data" is key for likelihood.

MCMC Inference Button

with coin_model:step = pm.Metropolis()      # focus on this, the Inference Button:coin_trace = pm.sample(2000, step=step)

结果：

pm.traceplot(coin_trace)
plt.show()

pm.plot_posterior(coin_trace[100:], color='#87ceeb',rope=[0.48,0.52], point_estimate='mean', ref_val=0.5)
plt.show()

例子2：化学活性问题

我有一个新开发的分子X； X在阻止流感方面的效果有多好？

实验

测试X的浓度范围，测量流感活动
根据实验结果计算 IC50：导致病毒复制率减半的X浓度。

数据

import numpy as np
import pandas as pdchem_data = [(0.00080, 99),
(0.00800, 91),
(0.08000, 89),
(0.40000, 89),
(0.80000, 79),
(1.60000, 61),
(4.00000, 39),
(8.00000, 25),
(80.00000, 4)]chem_df = pd.DataFrame(chem_data)
chem_df.columns = ['concentration', 'activity']
chem_df['concentration_log'] = chem_df['concentration'].apply(lambda x:np.log10(x))

参数问题

给出数据，化学品的IC50 值是多少, 以及其周围的不确定性？

说明

可视化数据

def plot_chemical_data(log=True):fig = plt.figure(figsize=(10,6))ax = fig.add_subplot(1,1,1)       if log:ax.scatter(x=chem_df['concentration_log'], y=chem_df['activity'])ax.set_xlabel('log10(concentration (mM))', fontsize=20)    else:ax.scatter(x=chem_df['concentration'], y=chem_df['activity'])ax.set_xlabel('concentration (mM)', fontsize=20)ax.set_xticklabels([int(i) for i in ax.get_xticks()], fontsize=18)ax.set_yticklabels([int(i) for i in ax.get_yticks()], fontsize=18)plt.hlines(y=50, xmin=min(ax.get_xlim()), xmax=max(ax.get_xlim()), linestyles='--',)    return figfig = plot_chemical_data(log=True)
plt.show()

with pm.Model() as ic50_model:beta = pm.HalfNormal('beta', sd=100**2)ic50_log10 = pm.Flat('IC50_log10')  # Flat prior# MATH WITH DISTRIBUTION OBJECTS!measurements = beta / (1 + np.exp(chem_df['concentration_log'].values - ic50_log10))y_like = pm.Normal('y_like', mu=measurements, observed=chem_df['activity'])  ic50 = pm.Deterministic('IC50', np.power(10, ic50_log10))# MCMC Inference Button
with ic50_model:step = pm.Metropolis()ic50_trace = pm.sample(10000, step=step)

pm.traceplot(ic50_trace[2000:], varnames=['IC50_log10', 'IC50'])  # live: sample from step 2000 onwards.
plt.show()

pm.plot_posterior(ic50_trace[4000:], varnames=['IC50'], color='#87ceeb', point_estimate='mean')
plt.show()

该化学物质的IC50在约 [2mM，2.4mM]（95％HPD）

问题类型2：实验组之间的比较

实验组和对照组的不同

例子1：药物IQ问题

药物治疗是否影响 IQ Scores

数据（包括对照实验数据）

drug = [  99.,  110.,  107.,  104., 省略]
placebo = [  95.,  105.,  103.,   99., 省略]def ECDF(data):x = np.sort(data)y = np.cumsum(x) / np.sum(x)        return x, ydef plot_drug():fig = plt.figure()ax = fig.add_subplot(1,1,1)x_drug, y_drug = ECDF(drug)ax.plot(x_drug, y_drug, label='drug, n={0}'.format(len(drug)))x_placebo, y_placebo = ECDF(placebo)ax.plot(x_placebo, y_placebo, label='placebo, n={0}'.format(len(placebo)))ax.legend()ax.set_xlabel('IQ Score')ax.set_ylabel('Cumulative Frequency')ax.hlines(0.5, ax.get_xlim()[0], ax.get_xlim()[1], linestyle='--')        return figfig = plot_drug()
plt.show()

from scipy.stats import ttest_indttest_ind(drug, placebo)# Ttest_indResult(statistic=2.2806701634329549, pvalue=0.025011500508647616)

实验

随机将参与者分配给两个实验组：
- +drug vs. -drug
测量每个参与者的 IQ Scores

说明

建模：

y_vals = np.concatenate([drug, placebo])
labels = ['drug'] * len(drug) + ['placebo'] * len(placebo)data = pd.DataFrame([y_vals, labels]).T
data.columns = ['IQ', 'treatment']
with pm.Model() as kruschke_model: # Linking Distribution Objects together is done by # passing objects into other objects' parameters.mu_drug = pm.Normal('mu_drug', mu=0, sd=100**2)mu_placebo = pm.Normal('mu_placebo', mu=0, sd=100**2)sigma_drug = pm.HalfCauchy('sigma_drug', beta=100)sigma_placebo = pm.HalfCauchy('sigma_placebo', beta=100)nu = pm.Exponential('nu', lam=1/29) + 1drug_like = pm.StudentT('drug', nu=nu, mu=mu_drug, sd=sigma_drug, observed=drug)placebo_like = pm.StudentT('placebo', nu=nu, mu=mu_placebo, sd=sigma_placebo, observed=placebo)diff_means = pm.Deterministic('diff_means', mu_drug - mu_placebo)pooled_sd = pm.Deterministic('pooled_sd', np.sqrt(np.power(sigma_drug, 2) + np.power(sigma_placebo, 2) / 2))effect_size = pm.Deterministic('effect_size', diff_means / pooled_sd)with kruschke_model:kruschke_trace = pm.sample(10000, step=pm.Metropolis())

结果：

pm.traceplot(kruschke_trace[2000:], varnames=['mu_drug', 'mu_placebo'])
plt.show()

pm.plot_posterior(kruschke_trace[2000:], color='#87ceeb',varnames=['mu_drug', 'mu_placebo', 'diff_means'])
plt.show()

用PyMC3进行贝叶斯统计分析（代码+实例）相关推荐

Adaboost算法原理分析和实例+代码（简明易懂）
Adaboost算法原理分析和实例+代码(简明易懂) [尊重原创,转载请注明出处] http://blog.csdn.net/guyuealian/article/details/70995333 ...
【Java 代码实例 14】BeanUtils用法详解，附源码分析
目录一.org.apache.commons.beanutils.BeanUtils简介二.使用的前置条件三.添加pom 四.org.apache.commons.beanutils.BeanU ...
使用python进行贝叶斯统计分析
本文讲解了使用PyMC3进行基本的贝叶斯统计分析过程. 最近我们被客户要求撰写关于贝叶斯统计分析的研究报告,包括一些图形和统计输出. # 导入 import pymc3 as pm # python的 ...
java从尾到头打印链表数据_Java编程实现从尾到头打印链表代码实例
问题描述:输入一个链表的头结点,从尾巴到头反过来打印出每个结点的值. 首先定义链表结点 public class ListNode { int val; ListNode next = null; L ...
mysql调试事件_mysql日志管理分析调试实例_mysql
以下的文章主要介绍的是mysql 操作日志查看的实际操作步骤以及对其实际操作步骤的具体描述,假如你在实际操作中遇到相似的情况,但是你却不知道对其如何正确的解决,那么以下的文章对你而言一定是良师益友. ...
ActiveMQ入门系列二：入门代码实例（点对点模式）
在上一篇<ActiveMQ入门系列一:认识并安装ActiveMQ(Windows下)>中,大致介绍了ActiveMQ和一些概念,并下载.安装.启动他,还访问了他的控制台页面. 这篇,就用代 ...
php面向对象代码_PHP面向对象之抽象类详解（代码实例）
[摘要] PHP即"超文本预处理器",是一种通用开源脚本语言.PHP是在服务器端执行的脚本语言,与C语言类似,是常用的网站编程语言.PHP独特的语法混合了C.Java.Perl以及 ...
编译原理中词法分析的递归下降分析法实例--能被5整除的二进制数---c语言实现
一.前言又到了一周一度的编译原理实验课,一次实验课上完了,又是大学生必备技能-写实验报告.行了,废话不多说,我直接展现,如何实现编译原理中词法分析的递归下降分析法实例–能被5整除的二进制数的思路.作 ...
python 监控linux硬盘,Python3监控windows,linux系统的CPU、硬盘、内存使用率和各个端口的开启情况详细代码实例...
由于项目的需要,需要做一个简单监控服务器的CPU利用率.CPU负载.硬盘使用率.内存利用率和服务器的各个端口的开启情况的程序,并把结果通知到监控平台,如果出现异常,监控平台打电话或者发短信通知给具体的 ...

用PyMC3进行贝叶斯统计分析（代码+实例）

例子1：抛硬币问题

参数问题

说明

例子2：化学活性问题

实验

数据

参数问题

说明

可视化数据

问题类型2：实验组之间的比较

例子1：药物IQ问题

实验

说明

用PyMC3进行贝叶斯统计分析（代码+实例）相关推荐

最新文章

热门文章