python - 使用 Pomegranate 拟合 Beta 分布

Question

我正在尝试使用库 pomegranate 来近似 Beta 分布。然而，当我尝试从生成的数据中近似参数时，我得到了非常不同的参数。重现该错误的代码如下

import numpy as np
from pomegranate import * 

X = np.random.beta(1, 5, size=10000).reshape(-1, 1) # sample from beta distribution with alpha = 1, beta = 5
print(BetaDistribution.from_samples(X).parameters) # approximate beta parameters
>>> [0.0, 10000.0] # error here

我不确定错误从何而来。好像是测试文件test_distributions.py产生正确的答案。如果有任何关于如何修复 pomegranate 或在 pomegranate 中创建自定义模型的建议，我们将不胜感激。

注意我使用的是Python 3.6.8

Answer 1

按照这个issue回答 , 当前库中提供的 BetaDistribution 是 beta 二项式分布而不是 beta 分布。这就是为什么该模型无法拟合 beta 分布的样本。

解决方案

我使用 BayesianOptimization 得到了解决方法图书馆。基本上，我尝试使用贝叶斯优化库最大化给定数据的分布的对数似然。以下代码对于混合分布也具有很好的概括性。

from bayes_opt import BayesianOptimization

data = np.random.beta(1, 5, size=10000) # create data

def beta_loss(a, b):
    beta_loss = BetaDistribution(a, b).probability(data)
    return np.log(beta_loss).sum()

optimizer = BayesianOptimization(
    f=beta_loss, 
    pbounds={'a': (0.5, 5), 
             'b': (0.5, 20)}, 
    random_state=10
)
# optimize the parameters
optimizer.maximize(
    init_points=5, 
    n_iter=100
)

# plot approximated distribution vs. distribution of the data
x = np.arange(0, 1, 0.01)
plt.hist(data, density=True, bins=100, alpha=0.1)
a, b = [v for k, v in optimizer.max['params'].items()]
plt.plot(x, BetaDistribution(a, b).probability(x))
plt.show()

额外(用于混合分布)

这里我只是举例说明如何优化Beta分布和高斯分布混合的参数:

from bayes_opt import BayesianOptimization

# example data of beta/gaussian distribution
data = np.hstack((np.random.beta(1, 10, size=2000), 
                  np.random.randn(1000) * 0.2 + 0.6))
data = data[np.logical_and(data >= 0.0, data <= 1.0)]

def loss_bimodal(a, b, mu, sigma, w1):
    beta_loss = BetaDistribution(a, b).probability(data)
    norm_loss = NormalDistribution(mu, sigma).probability(data)
    return np.log(w1 * beta_loss + (1 - w1) * norm_loss).sum()

def pdf_bimodal(a, b, mu, sigma, w1, x=np.arange(0, 1, 0.01)):
    return w1 * BetaDistribution(a, b).probability(x) + \
        (1 - w1) * NormalDistribution(mu, sigma).probability(x)

optimizer = BayesianOptimization(
    f=loss_bimodal, 
    pbounds={'mu': (0., 1.), 
             'sigma': (0., 1.), 
             'a': (0.5, 5), 
             'b': (1, 25), 
             'w1': (0., 1.)},
    random_state=1
)
optimizer.maximize(
    init_points=5, 
    n_iter=100
)

使用优化后的参数绘制分布图，如下所示:

a, b, mu, sigma, w1 = [v for k, v in optimizer.max['params'].items()]
x = np.arange(0, 1, 0.01)
plt.plot(x, pdf(a, b, mu, sigma, w1, x))
plt.hist(data, density=True, bins=100)
plt.show()