海外如何淘宝网站建设,2022新闻大事件摘抄,江苏省工程造价信息网,厦门市机场建设招投标网站问题类型1#xff1a;参数估计
真实值是否等于X#xff1f;
给出数据#xff0c;对于参数#xff0c;可能的值的概率分布是多少#xff1f;
例子1#xff1a;抛硬币问题
硬币扔了n次#xff0c;正面朝上是h次。
参数问题
想知道 p 的可能性。给定 n 扔的次数和 h … 问题类型1参数估计
真实值是否等于X
给出数据对于参数可能的值的概率分布是多少
例子1抛硬币问题
硬币扔了n次正面朝上是h次。
参数问题
想知道 p 的可能性。给定 n 扔的次数和 h 正面朝上次数p 的值很可能接近 0.5比如说在 [0.480.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, pp_prior, observedtosses) # observeddata is key for likelihood.
MCMC Inference Button
with coin_model:step pm.Metropolis() # focus on this, the Inference Button:coin_trace pm.sample(2000, stepstep)
结果
pm.traceplot(coin_trace)
plt.show() pm.plot_posterior(coin_trace[100:], color#87ceeb,rope[0.48,0.52], point_estimatemean, ref_val0.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(logTrue):fig plt.figure(figsize(10,6))ax fig.add_subplot(1,1,1) if log:ax.scatter(xchem_df[concentration_log], ychem_df[activity])ax.set_xlabel(log10(concentration (mM)), fontsize20) else:ax.scatter(xchem_df[concentration], ychem_df[activity])ax.set_xlabel(concentration (mM), fontsize20)ax.set_xticklabels([int(i) for i in ax.get_xticks()], fontsize18)ax.set_yticklabels([int(i) for i in ax.get_yticks()], fontsize18)plt.hlines(y50, xminmin(ax.get_xlim()), xmaxmax(ax.get_xlim()), linestyles--,) return figfig plot_chemical_data(logTrue)
plt.show() with pm.Model() as ic50_model:beta pm.HalfNormal(beta, sd100**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, mumeasurements, observedchem_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, stepstep)
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_estimatemean)
plt.show() 该化学物质的IC50在约 [2mM2.4mM]95HPD
问题类型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, labeldrug, n{0}.format(len(drug)))x_placebo, y_placebo ECDF(placebo)ax.plot(x_placebo, y_placebo, labelplacebo, 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(statistic2.2806701634329549, pvalue0.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, mu0, sd100**2)mu_placebo pm.Normal(mu_placebo, mu0, sd100**2)sigma_drug pm.HalfCauchy(sigma_drug, beta100)sigma_placebo pm.HalfCauchy(sigma_placebo, beta100)nu pm.Exponential(nu, lam1/29) 1drug_like pm.StudentT(drug, nunu, mumu_drug, sdsigma_drug, observeddrug)placebo_like pm.StudentT(placebo, nunu, mumu_placebo, sdsigma_placebo, observedplacebo)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, steppm.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()
