apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-06-26 15:49:40 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Use solve in lti_disc

Browse source
This commit is contained in:
Arno Solin 2013-11-11 22:53:45 +00:00
parent f4425efdae
commit 53f0fbcbff
1 changed files with 15 additions and 9 deletions
Download patch file Download diff file Expand all files Collapse all files

20
GPy/models/state_space.py
View file

@ -48,7 +48,7 @@ class StateSpace(Model):
self.ensure_default_constraints() self.ensure_default_constraints()
# Assert that the kernel is supported # Assert that the kernel is supported
#assert self.kern.sde(), "This kernel is not supported for state space estimation" #assert self.kern.sde() not False, "This kernel is not supported for state space estimation"
def _set_params(self, x): def _set_params(self, x):
self.kern._set_params(x[:self.kern.num_params_transformed()]) self.kern._set_params(x[:self.kern.num_params_transformed()])
@ -97,8 +97,6 @@ class StateSpace(Model):
# Run the Rauch-Tung-Striebel smoother # Run the Rauch-Tung-Striebel smoother
(M, P) = self.rts_smoother(F,L,Qc,X.T,M,P) (M, P) = self.rts_smoother(F,L,Qc,X.T,M,P)
#raise NameError('HiThere')
# Put the data back in the original order # Put the data back in the original order
M = M[:,return_inverse] M = M[:,return_inverse]
P = P[:,:,return_inverse] P = P[:,:,return_inverse]
@ -193,8 +191,15 @@ class StateSpace(Model):
return Y.T return Y.T
def posterior_samples(self, X, size=10): def posterior_samples(self, X, size=10):
# TODO
return self.posterior_samples_f(X,size) # Make samples of f
Y = self.posterior_samples_f(X,size)
# Add noise
Y += np.sqrt(self.sigma2)*np.random.randn(Y.shape)
# Return trajectory
return Y
def kalman_filter(self,F,L,Qc,H,R,Pinf,X,Y): def kalman_filter(self,F,L,Qc,H,R,Pinf,X,Y):
# KALMAN_FILTER - Run the Kalman filter for a given model and data # KALMAN_FILTER - Run the Kalman filter for a given model and data
@ -299,7 +304,7 @@ class StateSpace(Model):
f = np.zeros((F.shape[0],X.shape[1])) f = np.zeros((F.shape[0],X.shape[1]))
# Initial state # Initial state
f[:,0:1] = np.dot(np.linalg.cholesky(Pinf),np.random.randn(F.shape[0],1)) f[:,0:1] = np.linalg.cholesky(Pinf).dot(np.random.randn(F.shape[0],1))
# Sweep through remaining time points # Sweep through remaining time points
for k in range(1,X.shape[1]): for k in range(1,X.shape[1]):
@ -325,7 +330,8 @@ class StateSpace(Model):
Phi[:n,n:] = L.dot(Qc).dot(L.T) Phi[:n,n:] = L.dot(Qc).dot(L.T)
Phi[n:,n:] = -F.T Phi[n:,n:] = -F.T
AB = linalg.expm(Phi*dt).dot(np.vstack((np.zeros((n,n)),np.eye(n)))) AB = linalg.expm(Phi*dt).dot(np.vstack((np.zeros((n,n)),np.eye(n))))
Q = AB[:n,:].dot(linalg.inv(AB[n:,:])) #Q = AB[:n,:].dot(linalg.inv(AB[n:,:]))
Q = linalg.solve(AB[n:,:].T,AB[:n,:].T)
# The dynamical model # The dynamical model
A = linalg.expm(F*dt) A = linalg.expm(F*dt)