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

xt

Browse source
This commit is contained in:
mu 2013-11-29 17:14:49 +00:00
parent 71891a95d3
commit 1e01d79072
1 changed files with 35 additions and 18 deletions
Download patch file Download diff file Expand all files Collapse all files

53
GPy/models/state_space_xt_sep.py
View file

@ -45,6 +45,7 @@ class StateSpace_1(Model):
# Default kernel # Default kernel
if kernel is None: if kernel is None:
self.kern = kern.Matern32(1) self.kern = kern.Matern32(1)
self.spacekern = kern.Matern32(1)
else: else:
self.kern = kernel self.kern = kernel
@ -95,20 +96,23 @@ class StateSpace_1(Model):
# Get the model matrices from the kernel # Get the model matrices from the kernel
(F,L,Qc,H,Pinf) = self.kern.sde() (F,L,Qc,H,Pinf) = self.kern.sde()
#F = np.kron(eye(n),F1); n=X.shape[0]
#L = np.kron(eye(n),L1); F1 = np.kron(np.eye(n),F);
#Qc = K*Qc1; #kron(K,Qc1); L1 = np.kron(np.eye(n),L);
#H = np.kron(eye(n),H1); K1=self.spacekern.K(X)
#Pinf = np.kron(K,Pinf1); Qc1 = K1*Qc; #kron(K,Qc1);
H1 = np.kron(np.eye(n),H);
Pinf1 = np.kron(K1,Pinf);
# Run the Kalman filter # Run the Kalman filter
(M, P) = self.kalman_filter(F,L,Qc,H,self.sigma2,Pinf,X.T,Y.T) #(M, P) = self.kalman_filter(F,L,Qc,H,self.sigma2,Pinf,X.T,Y.T)
(M, P) = self.kalman_filter(F1,L1,Qc1,H1,self.sigma2,Pinf1,X.T,Y.T)
# Run the Rauch-Tung-Striebel smoother # Run the Rauch-Tung-Striebel smoother
#if not filter: #if not filter:
(M, P) = self.rts_smoother(F,L,Qc,X.T,M,P) #(M, P) = self.rts_smoother(F,L,Qc,X.T,M,P)
#(M, P) = self.rts_smoother(F1,L1,Qc1,X.T,M,P)
# Put the data back in the original order # Put the data back in the original order
M = M[:,return_inverse] M = M[:,return_inverse]
@ -123,7 +127,8 @@ class StateSpace_1(Model):
V = np.tensordot(H[0],P,(0,0)) V = np.tensordot(H[0],P,(0,0))
V = np.tensordot(V,H[0],(0,0)) V = np.tensordot(V,H[0],(0,0))
V = V[:,None] V = V[:,None]
sfasdfasdfasdf=pppp
# Return the posterior of the state # Return the posterior of the state
return (m, V) return (m, V)
@ -154,6 +159,8 @@ class StateSpace_1(Model):
# Define the frame on which to plot # Define the frame on which to plot
resolution = resolution or 200 resolution = resolution or 200
Xgrid, xmin, xmax = x_frame1D(self.X, plot_limits=plot_limits) Xgrid, xmin, xmax = x_frame1D(self.X, plot_limits=plot_limits)
# T grid???
# Make a prediction on the frame and plot it # Make a prediction on the frame and plot it
if plot_raw: if plot_raw:
@ -244,17 +251,27 @@ class StateSpace_1(Model):
# Prediction step # Prediction step
MF[:,k] = A.dot(MF[:,k-1]) MF[:,k] = A.dot(MF[:,k-1])
PF[:,:,k] = A.dot(PF[:,:,k-1]).dot(A.T) + Q PF[:,:,k] = A.dot(PF[:,:,k-1]).dot(A.T) + Q
# Update step (only if there is data) # Update step (only if there is data)
#if not np.isnan(Y[:,k]):
# if Y.shape[0]==1:
# K = PF[:,:,k].dot(H.T)/(H.dot(PF[:,:,k]).dot(H.T) + R)
# else:
# LL = linalg.cho_factor(H.dot(PF[:,:,k]).dot(H.T) + R)
# K = linalg.cho_solve(LL, H.dot(PF[:,:,k].T)).T
# stop
# MF[:,k] += K.dot(Y[:,k]-H.dot(MF[:,k]))
# PF[:,:,k] -= K.dot(H).dot(PF[:,:,k])
if not np.isnan(Y[:,k]): if not np.isnan(Y[:,k]):
if Y.shape[0]==1:
K = PF[:,:,k].dot(H.T)/(H.dot(PF[:,:,k]).dot(H.T) + R) LL = linalg.cho_factor(H.dot(PF[:,:,k]).dot(H.T) + R)
else: K = linalg.cho_solve(LL, H.dot(PF[:,:,k].T)).T
LL = linalg.cho_factor(H.dot(PF[:,:,k]).dot(H.T) + R)
K = linalg.cho_solve(LL, H.dot(PF[:,:,k].T)).T
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])
stop
# Return values # Return values
return (MF, PF) return (MF, PF)