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

FIX: Fixe bug with "expm" function in "state_space_new". Also some minor changes

Browse source 
Test function has been modified also.
This commit is contained in:
Alexander Grigorievskiy 2015-04-08 19:44:19 +03:00
parent 5ff256079b
commit d9cf9c3bff
2 changed files with 29 additions and 22 deletions
Download patch file Download diff file Expand all files Collapse all files

38
GPy/examples/state_space.py
View file

@ -7,36 +7,36 @@ import GPy.models.state_space_new as SS_new
#X = np.linspace(0, 10, 2000)[:, None] #X = np.linspace(0, 10, 2000)[:, None]
#Y = np.sin(X) + np.random.randn(*X.shape)*0.1 #Y = np.sin(X) + np.random.randn(*X.shape)*0.1
## Need to run these lines when X and Y are imported -> # Need to run these lines when X and Y are imported ->
#X.shape = (X.shape[0],1) X.shape = (X.shape[0],1)
#Y.shape = (Y.shape[0],1) Y.shape = (Y.shape[0],1)
## Need to run these lines when X and Y are imported <- # Need to run these lines when X and Y are imported <-
# Generation of minimal example data -> ## Generation of minimal example data ->
X = np.random.rand(3) #X = np.random.rand(3)
sort_index = np.argsort(X) #sort_index = np.argsort(X)
X = X[sort_index]; X.shape = (X.shape[0],1) #X = X[sort_index]; X.shape = (X.shape[0],1)
Y = np.sin(10*X) + np.random.randn(*X.shape)*0.1 #Y = np.sin(10*X) + np.random.randn(*X.shape)*0.1
# Generation of minimal example data <- ## Generation of minimal example data <-
#plt.figure() #plt.figure()
#plt.plot( X, Y) #plt.plot( X, Y)
#plt.show() #plt.show()
kernel = GPy.kern.Matern32(X.shape[1]) #kernel = GPy.kern.Matern32(X.shape[1])
m = GPy.models.StateSpace(X,Y, kernel) #m = GPy.models.StateSpace(X,Y, kernel)
print m
# #
m.optimize() #print m
# ##
print m #m.optimize(optimizer='bfgs',messages=True)
##
#print m
kernel1 = GPy.kern.Matern32(X.shape[1]) kernel1 = GPy.kern.Matern32(X.shape[1])
m1 = GPy.models.GPRegression(X,Y, kernel1) m1 = GPy.models.GPRegression(X,Y, kernel1)
print m1 print m1
m1.optimize() m1.optimize(optimizer='bfgs',messages=True)
print m1 print m1
@ -45,7 +45,7 @@ m2 = SS_new.StateSpace(X,Y, kernel2)
print m2 print m2
m2.optimize() m2.optimize(optimizer='bfgs',messages=True)
print m2 print m2

13
GPy/models/state_space_main.py
View file

@ -657,6 +657,7 @@ class DescreteStateSpace(object):
dP_pred[:,:,j] += dP_pred[:,:,j].T dP_pred[:,:,j] += dP_pred[:,:,j].T
dP_pred[:,:,j] += np.dot( A ,np.dot(dP, A.T)) + dQ dP_pred[:,:,j] += np.dot( A ,np.dot(dP, A.T)) + dQ
dP_pred[:,:,j] = 0.5*(dP_pred[:,:,j] + dP_pred[:,:,j].T) #symmetrize
else: else:
dm_pred = None dm_pred = None
dP_pred = None dP_pred = None
@ -1439,7 +1440,7 @@ class ContDescrStateSpace(DescreteStateSpace):
unique_round_decimals = 8 unique_round_decimals = 8
dt = np.empty((X.shape[0],)) dt = np.empty((X.shape[0],))
dt[1:] = np.diff(X[:,0],axis=0) dt[1:] = np.diff(X[:,0],axis=0)
dt[0] = dt[1] dt[0] = 0#dt[1]
unique_indices = np.unique(np.round(dt, decimals=unique_round_decimals)) unique_indices = np.unique(np.round(dt, decimals=unique_round_decimals))
if len(unique_indices) > 20: if len(unique_indices) > 20:
@ -1531,13 +1532,19 @@ class ContDescrStateSpace(DescreteStateSpace):
# The dynamical model # The dynamical model
A = linalg.expm(F*dt) A = linalg.expm(F*dt)
if np.any( np.isnan(A)):
A = linalg.expm3(F*dt)
# The covariance matrix Q by matrix fraction decomposition -> # The covariance matrix Q by matrix fraction decomposition ->
Phi = np.zeros((2*n,2*n)) Phi = np.zeros((2*n,2*n))
Phi[:n,:n] = F Phi[:n,:n] = F
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)
if np.any(np.isnan(AB)):
AB = linalg.expm3(Phi*dt)
AB = np.dot(AB, np.vstack((np.zeros((n,n)),np.eye(n))))
Q_noise_1 = linalg.solve(AB[n:,:].T,AB[:n,:].T) Q_noise_1 = linalg.solve(AB[n:,:].T,AB[:n,:].T)
# The covariance matrix Q by matrix fraction decomposition <- # The covariance matrix Q by matrix fraction decomposition <-