【发布时间】:2014-05-22 09:37:27
【问题描述】:
我正在 Python/NumPy 中实现贝叶斯变更点检测(如果您有兴趣,请查看 paper)。我需要计算[a, b] 范围内数据的可能性,其中a 和b 可以具有从1 到n 的所有值。但是,我可以在某些点修剪计算,这样我就不必计算所有可能性。另一方面,一些可能性被多次使用,因此我可以通过将值保存在矩阵P[a, b] 中来节省时间。现在,每当我使用它时,我都会检查该值是否已经计算过,但我发现这有点麻烦。它看起来像这样:
# ...
P = np.ones((n, n)) * np.inf # a likelihood can't get inf, so I use it
# as pseudo value
for a in range(n):
for b in range(a, n):
# The following two lines get annoying and error prone if you
# use P more than once
if P[a, b] == np.inf:
P[a, b] = likelihood(data, a, b)
Q[a] += P[a, b] * g[a] * Q[a - 1] # some computation using P[a, b]
# ...
我想知道,是否有更直观和 Pythonic 的方式来实现这一点,而无需在每次使用 P[a, b] 之前使用 if ... 语句。如果不满足某些条件,则类似于自动函数调用。我当然可以让likelihood 函数意识到它可以保存值的事实,但是它需要某种状态(例如成为一个对象)。我想避免这种情况。
似然函数
由于在评论中要求,我添加了似然函数。它实际上先计算共轭,然后计算可能性。而且都是对数表示...所以它相当复杂。
from scipy.special import gammaln
def gaussian_obs_log_likelihood(data, t, s):
n = s - t
mean = data[t:s].sum() / n
muT = (n * mean) / (1 + n)
nuT = 1 + n
alphaT = 1 + n / 2
betaT = 1 + 0.5 * ((data[t:s] - mean) ** 2).sum() + ((n)/(1 + n)) * (mean**2 / 2)
scale = (betaT*(nuT + 1))/(alphaT * nuT)
# splitting the PDF of the student distribution up is /much/ faster. (~ factor 20)
prob = 1
for yi in data[t:s]:
prob += np.log(1 + (yi - muT)**2/(nuT * scale))
lgA = gammaln((nuT + 1) / 2) - np.log(np.sqrt(np.pi * nuT * scale)) - gammaln(nuT/2)
return n * lgA - (nuT + 1)/2 * prob
虽然我使用 Python 2.7,但对 2.7 和 3.x 的答案都表示赞赏。
【问题讨论】:
-
你能告诉我们你的似然函数吗?使该函数适用于向量而不是标量
-
可能是初始值
None?然后你可以使用P[a, b] = likelihood(data, a, b) if not P[a, b] -
它实际上适用于向量。
a和b只是计算可能性的data的范围。反正我会加的。 -
既然整个数组
P被初始化为inf,那么if的条件不应该总是为真吗?或者您是否在代码的其他部分更改P? -
P[a, b]被多次使用。所以if...语句只用于计算一次。第一次需要P[a, b]。我实际上在if语句中更改了它。代码做了一些简化,让问题更清晰。
标签: python numpy conditional-statements dynamic-programming