贝叶法则 - 如何计算可能性



给定一些数据data,它对应于一个二元掷硬币序列,其中正面为1,反面为0。θ是一个介于0到1之间的值,表示硬币被抛掷时正面朝上的概率。

如何计算可能性?我依稀记得一个公式:

likelihood = (theta)^(h)*(1-theta)^(1-h) 

,其中正面h = 1,反面h = 0。我实现了以下代码:

import numpy as np
(np.prod([theta*1 for i in data if i==1]) * np.prod([1-theta for i in data if i==0]))

这段代码适用于某些情况,但不适用于一些隐藏的情况(所以我不确定它有什么问题)。

有两种方法可以解释您正在尝试计算的内容:

  1. 该序列的概率,包括头部出现的顺序(这就是你在这里提出问题的方式)
  2. 在你的序列中出现的正面数量(我们称之为X)的概率,无论顺序如何(这是我认为你想要的)。

选项1:

import numpy as np
theta = 0.2      # Probability of H is 0.2, hence NOT a fair coin
data = [0, 1, 0, 1, 1, 1, 0, 0, 1, 1]     # T, H, T, H, H, .... 

def likelihood(theta, h):
    return (theta)**(h)*(1-theta)**(1-h)
likelihood(theta, 1) # 0.2
likelihood(theta, 0) # 0.8
singlethrow = [likelihood(theta, x) for x in data]
prob1 = np.prod(singlethrow)    # 2.6214400000000015e-05

prob1将很快收敛于零,因为每次额外的抛硬币都会将现有的概率乘以小于1的数字(正面为0.2,反面为0.8)

选项2:

是二项分布。这是所有可能结果的概率相加,比如抛10次硬币得到6次正面。我们已经在上面的选项1中评估了10次投掷6次正面的特定序列。有210种这样的方法(= 10!/(6!*(10−6)!))

scipy.stats.binom.pmf()功能为您计算此概率:

import scipy, scipy.stats
prob2 = scipy.stats.binom.pmf(6, 10, theta)

或者,更一般地,如果您依赖于我上面定义的形式的data:

X = sum([toss == 1 for toss in data])
N = len(data)
prob3 = scipy.stats.binom.pmf(X, N, theta)
prob2 == prob3   # True 

如果您对贝叶斯方法感兴趣,您可能想看看conjugate_prior

from conjugate_prior import BetaBinomial
prior_model = BetaBinomial(1,1) # Uninformative prior
updated_model = prior_model.update(heads, tails)
credible_interval = updated_model.posterior(0.45, 0.55)
print ("There's {p:.2f}% chance that the coin is fair".format(p=credible_interval*100))

最新更新