根据伯努利分布的定义，我们可以得到以下公式： f ( x ) = p x ( 1 − p ) 1 − x for  x ∈ { 0 , 1 } f(x)=p^{x}(1-p)^{1-x} \quad \text { for } x \in\{0,1\} f(x)=px(1−p)1−x for x∈{0,1}当我们进行观测，得到了事件$x$发生频率的观测值（代码中假设是0.5，比如抛硬币1000次，500次朝上、500次朝下），我们就可以画出伯努利分布的概率p的密度图。

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
#x表示观测到的实验结果，也就是有0.5的概率正面朝上
x = 0.5
p = np.linspace(0,1,100)
f = (p**x)*((1-p)**(1-x))#以下使用spline对曲线进行平滑处理
import os
os.sys.path.append('p:\\users\\loara\anaconda3\\lib\\site-packages')
from scipy.interpolate import make_interp_spline
p_smooth = np.linspace(p.min(),p.max(),300) #300 represents number of points to make between T.min and T.max
f_smooth = make_interp_spline(p,f)(p_smooth)
plt.plot(p_smooth,f_smooth)
plt.show()

可见，当$p$取0.5时，概率密度达到最大，所以我们认为对$p$的最佳估计就是0.5。

我们来修改一下实验数据，如果观测到的$x$为0.7，我们可以得到以下图形：

我们可以看到最佳估计就在0.7估计，此时此刻，我只能说：伯努利分布真牛逼！