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

sep xt

Browse source
This commit is contained in:
mu 2013-12-09 17:59:33 +00:00
parent 94c7047b82
commit d44e16e615
1 changed files with 63 additions and 41 deletions
Download patch file Download diff file Expand all files Collapse all files

84
GPy/models/state_space_xt_sep.py
View file

@ -22,22 +22,22 @@ from GPy.util.plot import gpplot, Tango, x_frame1D
import pylab as pb import pylab as pb
class StateSpace_1(Model): class StateSpace_1(Model):
def __init__(self, X, Y, kernel=None): def __init__(self, Sp,X, Y, kernel=None):
super(StateSpace_1, self).__init__() super(StateSpace_1, self).__init__()
self.num_data, input_dim = X.shape self.num_data, input_dim = X.shape
assert input_dim==1, "State space methods for time and space 2" assert input_dim==1, "State space methods for time and space 2"
num_data_Y, self.output_dim = Y.shape num_data_Y, self.output_dim = Y.shape
assert num_data_Y == self.num_data, "X and Y data don't match" assert num_data_Y == self.num_data, "X and Y data don't match"
assert self.output_dim == 1, "State space methods for single outputs only" #assert self.output_dim == 1, "State space methods for single outputs only"
# Make sure the observations are ordered in time # Make sure the observations are ordered in time
sort_index = np.argsort(X[:,0]) sort_index = np.argsort(X[:,0])
self.X = X[sort_index] self.X = X[sort_index]
self.Y = Y[sort_index] self.Y = Y[sort_index]
sort_index = np.argsort(X[:,0]) #sort_index = np.argsort(X[:,0])
self.X = X[sort_index] #self.X = X[sort_index]
self.Y = Y[sort_index] #self.Y = Y[sort_index]
# Noise variance # Noise variance
self.sigma2 = 1. self.sigma2 = 1.
@ -45,7 +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) self.spacekern = kern.rbf(1)
else: else:
self.kern = kernel self.kern = kernel
@ -85,8 +85,10 @@ class StateSpace_1(Model):
def predict_raw(self, Xnew, filteronly=False): def predict_raw(self, Xnew, filteronly=False):
# Make a single matrix containing training and testing points # Make a single matrix containing training and testing points
X = np.vstack((self.X, Xnew)) #X = np.vstack((self.X, Xnew))
Y = np.vstack((self.Y, np.nan*np.zeros(Xnew.shape))) #Y = np.vstack((self.Y, np.nan*np.zeros(Xnew.shape)))
X=self.X
Y=self.Y
# Sort the matrix (save the order) # Sort the matrix (save the order)
_, return_index, return_inverse = np.unique(X,True,True) _, return_index, return_inverse = np.unique(X,True,True)
@ -97,17 +99,17 @@ class StateSpace_1(Model):
(F,L,Qc,H,Pinf) = self.kern.sde() (F,L,Qc,H,Pinf) = self.kern.sde()
n=X.shape[0] n=X.shape[0]
F1 = np.kron(np.eye(n),F); F1 = np.kron(np.eye(n),F)
L1 = np.kron(np.eye(n),L); L1 = np.kron(np.eye(n),L)
K1=self.spacekern.K(X) K1=self.spacekern.K(X)
Qc1 = K1*Qc; #kron(K,Qc1); Qc1 = K1*Qc #kron(K,Qc1);
H1 = np.kron(np.eye(n),H); H1 = np.kron(np.eye(n),H)
Pinf1 = np.kron(K1,Pinf); 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) (M, P) = self.kalman_filter(F1,L1,Qc1,H1,self.sigma2,Pinf1,X.T,Y)
# Run the Rauch-Tung-Striebel smoother # Run the Rauch-Tung-Striebel smoother
#if not filter: #if not filter:
@ -119,16 +121,20 @@ class StateSpace_1(Model):
P = P[:,:,return_inverse] P = P[:,:,return_inverse]
# Only return the values for Xnew # Only return the values for Xnew
M = M[:,self.num_data:] #M = M[:,self.num_data:]
P = P[:,:,self.num_data:] #P = P[:,:,self.num_data:]
# Calculate the mean and variance # Calculate the mean and variance
m = H.dot(M).T #m = H.dot(M).T
V = np.tensordot(H[0],P,(0,0)) m = H1.dot(M)
V = np.tensordot(V,H[0],(0,0)) #V1 = np.tensordot(H[0],P,(0,0))
V = V[:,None] #V2 = np.tensordot(V1,H[0],(0,0))
V=P[::F.shape[0],::F.shape[0],:]
#V1 = np.tensordot(H1.T,P,(0,0))
#V2 = np.tensordot(V1,H1,(1,1))
#stop
#V3 = V2[:,None]
sfasdfasdfasdf=pppp
# Return the posterior of the state # Return the posterior of the state
return (m, V) return (m, V)
@ -138,8 +144,8 @@ class StateSpace_1(Model):
(m, V) = self.predict_raw(Xnew,filteronly=filteronly) (m, V) = self.predict_raw(Xnew,filteronly=filteronly)
# Add the noise variance to the state variance # Add the noise variance to the state variance
V += self.sigma2 V += self.sigma2*np.eye(m.shape[0])
stop
# Lower and upper bounds # Lower and upper bounds
lower = m - 2*np.sqrt(V) lower = m - 2*np.sqrt(V)
upper = m + 2*np.sqrt(V) upper = m + 2*np.sqrt(V)
@ -164,16 +170,33 @@ class StateSpace_1(Model):
# Make a prediction on the frame and plot it # Make a prediction on the frame and plot it
if plot_raw: if plot_raw:
#m, v = self.predict_raw(Xgrid,filteronly=plot_filter)
m, v = self.predict_raw(Xgrid,filteronly=plot_filter) m, v = self.predict_raw(Xgrid,filteronly=plot_filter)
lower = m - 2*np.sqrt(v)
upper = m + 2*np.sqrt(v)
Y = self.Y Y = self.Y
#np.random.multivariant_normal(m[:,0],v[:,:,0])
#np.random.multivariant_normal(m[:,1],v[:,:,1])
pb.figure(1)
pb.imshow(m,interpolation="nearest")
pb.figure(2)
pb.imshow(Y,interpolation="nearest")
stop
#lower = m - 2*np.sqrt(v)
#upper = m + 2*np.sqrt(v)
else: else:
m, v, lower, upper = self.predict(Xgrid,filteronly=plot_filter) m, v, lower, upper = self.predict(Xgrid,filteronly=plot_filter)
Y = self.Y Y = self.Y
# Plot the values # Plot the values
gpplot(Xgrid, m, lower, upper, axes=ax, edgecol=linecol, fillcol=fillcol) #gpplot(Xgrid, m, lower, upper, axes=ax, edgecol=linecol, fillcol=fillcol)
gpplot(self.X, m, lower, upper, axes=ax, edgecol=linecol, fillcol=fillcol)
ax.plot(self.X, self.Y, 'kx', mew=1.5) ax.plot(self.X, self.Y, 'kx', mew=1.5)
# Optionally plot some samples # Optionally plot some samples
@ -247,7 +270,6 @@ class StateSpace_1(Model):
#(A, Q) = self.lti_disc(F,L,Qc,dt[:,k]) #(A, Q) = self.lti_disc(F,L,Qc,dt[:,k])
A = As[:,:,index[k]]; A = As[:,:,index[k]];
Q = Qs[:,:,index[k]]; Q = Qs[:,:,index[k]];
# 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
@ -263,14 +285,14 @@ class StateSpace_1(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])
if not np.isnan(Y[:,k]): #if not np.isnan(Y[:,k]):
LL = linalg.cho_factor(H.dot(PF[:,:,k]).dot(H.T) + R) LL = linalg.cho_factor(H.dot(PF[:,:,k]).dot(H.T) + R*np.eye(Y.shape[1]))
K = linalg.cho_solve(LL, H.dot(PF[:,:,k].T)).T K = linalg.cho_solve(LL, H.dot(PF[:,:,k].T)).T
stop
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)