给定一些数据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]))
这段代码适用于某些情况,但不适用于一些隐藏的情况(所以我不确定它有什么问题)。
有两种方法可以解释您正在尝试计算的内容:
- 该序列的概率,包括头部出现的顺序(这就是你在这里提出问题的方式)
- 在你的序列中出现的正面数量(我们称之为
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))