Fixed bug in log_likelihood evaluation

This commit is contained in:
Arno Solin 2013-11-10 20:38:41 +00:00
parent 8008df39eb
commit d6bae3364e

View file

@ -182,6 +182,8 @@ class StateSpace(Model):
MF[:,k] += K.dot(Y[:,k]-H.dot(MF[:,k])) MF[:,k] += K.dot(Y[:,k]-H.dot(MF[:,k]))
PF[:,:,k] -= K.dot(H).dot(PF[:,:,k]) PF[:,:,k] -= K.dot(H).dot(PF[:,:,k])
print K
# Return values # Return values
return (MF, PF) return (MF, PF)
@ -220,7 +222,7 @@ class StateSpace(Model):
dt = np.empty(X.shape) dt = np.empty(X.shape)
dt[:,0] = X[:,1]-X[:,0] dt[:,0] = X[:,1]-X[:,0]
dt[:,1:] = np.diff(X) dt[:,1:] = np.diff(X)
# Kalman filter for likelihood evaluation # Kalman filter for likelihood evaluation
for k in range(0,Y.shape[1]): for k in range(0,Y.shape[1]):
@ -235,14 +237,15 @@ class StateSpace(Model):
if not np.isnan(Y[:,k]): if not np.isnan(Y[:,k]):
v = Y[:,k]-H.dot(m) v = Y[:,k]-H.dot(m)
LL, isupper = linalg.cho_factor(H.dot(P).dot(H.T) + R) LL, isupper = linalg.cho_factor(H.dot(P).dot(H.T) + R)
lik -= 0.5*np.sum(np.log(np.diag(LL))) lik -= np.sum(np.log(np.diag(LL)))
lik -= 0.5*v.shape[0]*np.log(2*np.pi)
lik -= 0.5*linalg.cho_solve((LL, isupper),v).dot(v) lik -= 0.5*linalg.cho_solve((LL, isupper),v).dot(v)
K = linalg.cho_solve((LL, isupper), H.dot(P.T)).T K = linalg.cho_solve((LL, isupper), H.dot(P.T)).T
m += K.dot(v) m += K.dot(v)
P -= K.dot(H).dot(P) P -= K.dot(H).dot(P)
# Return likelihood # Return likelihood
return lik return lik[0,0]
def simulate(self,F,L,Qc,Pinf,X): def simulate(self,F,L,Qc,Pinf,X):
# Simulate a trajectory using the state space model # Simulate a trajectory using the state space model