Python KalmanFilter.filter_update方法代码示例

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
class StateEstimationAltitudeCam:
    def __init__(self):
        self.state = [0, 0]
        self.covariance = np.eye(2)
        self.observation_covariance = np.array([[1]])
        self.transition_covariance = np.array([[0.001, 0.001], [0.001, 0.001]])
        self.kf = KalmanFilter(
            transition_covariance=self.transition_covariance,  # H
            observation_covariance=self.observation_covariance,  # Q
        )
        self.previous_update = None

    def update(self, observations):
        if not self.previous_update:
            self.previous_update = time.time()
        dt = time.time() - self.previous_update

        self.state, self.covariance = self.kf.filter_update(
            self.state,
            self.covariance,
            observations,
            transition_matrix=np.array([[1, dt], [0, 1]]),
            observation_matrix=np.array([[1, 0]]),
            # observation_offset = np.array([, 0, 0])
            # observation_covariance=np.array(0.1*np.eye(1))
        )
        self.previous_update = time.time()

    def getAltitude(self):
        return self.state[0]

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
class KalmanMovingAverage(object):
    """
    Estimates the moving average of a price process via Kalman Filtering, using pykalman
    """

    def __init__(self, asset, observation_covariance=1.0, initial_value=0,
                 initial_state_covariance=1.0, transition_covariance=0.05,
                 initial_window=30, maxlen=300, freq='1d'):
        self.asset = asset
        self.freq = freq
        self.initial_window = initial_window

        self.kf = KalmanFilter(transition_matrices=[1],
                               observation_matrices=[1],
                               initial_state_mean=initial_value,
                               initial_state_covariance=initial_state_covariance,
                               observation_covariance=observation_covariance,
                               transition_covariance=transition_covariance)
        self.state_means = pd.Series([self.kf.initial_state_mean], name=self.asset)
        self.state_covs = pd.Series([self.kf.initial_state_covariance], name=self.asset)

    def update(self, observations):
        for dt, observation in observations[self.asset].iterkv():
            self._update(dt, observation)

    def _update(self, dt, observation):
        mu, cov = self.kf.filter_update(self.state_means.iloc[-1],
                                        self.state_covs.iloc[-1],
                                        observation)
        self.state_means[dt] = mu.flatten()[0]
        self.state_covs[dt] = cov.flatten()[0]

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
def test_kalman_filter_update():
    kf = KalmanFilter(
        data.transition_matrix,
        data.observation_matrix,
        data.transition_covariance,
        data.observation_covariance,
        data.transition_offsets,
        data.observation_offset,
        data.initial_state_mean,
        data.initial_state_covariance)

    # use Kalman Filter
    (x_filt, V_filt) = kf.filter(X=data.observations)

    # use online Kalman Filter
    n_timesteps = data.observations.shape[0]
    n_dim_state = data.transition_matrix.shape[0]
    x_filt2 = np.zeros((n_timesteps, n_dim_state))
    V_filt2 = np.zeros((n_timesteps, n_dim_state, n_dim_state))
    for t in range(n_timesteps - 1):
        if t == 0:
            x_filt2[0] = data.initial_state_mean
            V_filt2[0] = data.initial_state_covariance
        (x_filt2[t + 1], V_filt2[t + 1]) = kf.filter_update(
            x_filt2[t], V_filt2[t], data.observations[t + 1],
            transition_offset=data.transition_offsets[t]
        )
    assert_array_almost_equal(x_filt, x_filt2)
    assert_array_almost_equal(V_filt, V_filt2)

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
class KalmanRegression(object):
    """
    Uses a Kalman Filter to estimate regression parameters 
    in an online fashion.
    
    Estimated model: y ~ beta * x + alpha
    """
    
    def __init__(self, initial_y, initial_x, delta=1e-5, maxlen=3000):
        self._x = initial_x.name
        self._y = initial_y.name
        self.maxlen = maxlen
        trans_cov = delta / (1 - delta) * np.eye(2)
        obs_mat = np.expand_dims(
            np.vstack([[initial_x], [np.ones(initial_x.shape[0])]]).T, axis=1)
        
        self.kf = KalmanFilter(n_dim_obs=1, n_dim_state=2,
                               initial_state_mean=np.zeros(2),
                               initial_state_covariance=np.ones((2, 2)),
                               transition_matrices=np.eye(2),
                               observation_matrices=obs_mat,
                               observation_covariance=1.0,
                               transition_covariance=trans_cov)
        state_means, state_covs = self.kf.filter(initial_y.values)
        self.means = pd.DataFrame(state_means, 
                                  index=initial_y.index, 
                                  columns=['beta', 'alpha'])
        self.state_cov = state_covs[-1]
        
    def update(self, observations):
        x = observations[self._x]
        y = observations[self._y]
        mu, self.state_cov = self.kf.filter_update(
            self.state_mean, self.state_cov, y, 
            observation_matrix=np.array([[x, 1.0]]))
        mu = pd.Series(mu, index=['beta', 'alpha'], 
                       name=observations.name)
        self.means = self.means.append(mu)
        if self.means.shape[0] > self.maxlen:
            self.means = self.means.iloc[-self.maxlen:]
        
    def get_spread(self, observations):
        x = observations[self._x]
        y = observations[self._y]
        return y - (self.means.beta[-1] * x + self.means.alpha[-1])
        
    @property
    def state_mean(self):
        return self.means.iloc[-1]
        

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
def trackKalman():

    l_publ = rospy.Publisher("tracking", tracking, queue_size = 10)

    rate = rospy.Rate(10) # 10hz

    initstate = [current_measurement[0], current_measurement[1], 0, 0]
    Transition_Matrix=[[1,0,1,0],[0,1,0,1],[0,0,1,0],[0,0,0,1]]
    Observation_Matrix=[[1,0,0,0],[0,1,0,0]]
    initcovariance=1.0e-3*np.eye(4)
    transistionCov=1.0e-4*np.eye(4)
    observationCov=1.0e-1*np.eye(2)
    kf=KalmanFilter(transition_matrices=Transition_Matrix,
                observation_matrices =Observation_Matrix,
                initial_state_mean=initstate,
                initial_state_covariance=initcovariance,
                transition_covariance=transistionCov,
                observation_covariance=observationCov)

    start = 1
    t = 1
    while not rospy.is_shutdown():

        #cv2.rectangle(image,(face[0]-50,face[1]-50),(face[0]+50,face[1]+50),(255,0,0),2)
        ##cv2.rectangle(image,(face_center_pixels[0]-50,face_center_pixels[1]-50),(face_center_pixels[0]+50,face_center_pixels[1]+50),(255,0,0),2)
        #cv2.imshow("Calibrated Boundary", image)
        #cv2.waitKey(1)

        if (start == 1):
            start = 0
            filtered_state_means = initstate
            filtered_state_covariances = initcovariance
        
        print 'current measurement: ', current_measurement
        (pred_filtered_state_means, pred_filtered_state_covariances) = kf.filter_update(filtered_state_means, filtered_state_covariances, current_measurement);
        t += 1
        filtered_state_means = pred_filtered_state_means
        filtered_state_covariances = pred_filtered_state_covariances
        print 'predicted: ', filtered_state_means[0], filtered_state_means[1]
        #print type(current_measurement[0])
        #print type(filtered_state_means[0])
        location = tracking()
        location.x0 = current_measurement[0]
        location.y0 = current_measurement[1]
        location.x1 = filtered_state_means[0]
        location.y1 = filtered_state_means[1]
        l_publ.publish(location)
        print '\n'
        rate.sleep()

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
class Regression(object):
    """
    y = beta * x + alpha
    """

    def __init__(self, initial_x, initial_y, date, delta=1e-5):
        trans_cov = delta / (1 - delta) * np.eye(2)
        xList = []
        for x in initial_x:
            xList.append([[x,1.]])
        obs_mat = np.vstack([xList])

        self.kf = KalmanFilter(n_dim_obs=1, n_dim_state=2,
                               initial_state_mean=np.zeros(2),
                               initial_state_covariance=np.ones((2, 2)),
                               transition_matrices=np.eye(2),
                               observation_matrices=obs_mat,
                               observation_covariance=1.0,
                               transition_covariance=trans_cov)
        state_means, state_covs = self.kf.filter(np.array(initial_y))
        self.means = pd.DataFrame(state_means,
                                  columns=['beta', 'alpha'])
        self.state_cov = state_covs[-1]

    def update(self, observations, date):
        x = observations[0]
        y = observations[1]
        mu, self.state_cov = self.kf.filter_update(self.state_mean, self.state_cov, y,
                                                   observation_matrix=np.array([[x, 1.0]]))
        mu = pd.Series(mu, index=['beta', 'alpha'],
                       name=date)
        self.means = self.means.append(mu)

    def get_spread(self, observations):
        x = observations[0]
        y = observations[1]
        return y - (self.means.beta * x + self.means.alpha)

    @property
    def state_mean(self):
        return self.means.iloc[-1]

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
def kalman_filtering(price_sequence):

    h = 1 #time step
    A = np.array([[1,h,.5*h**2],
                [0,1,h],
                [0,0,1]])
                
    Q = np.eye(A.shape[0])
    #2) Apply the filter    
    kf = KalmanFilter(transition_matrices = A , transition_covariance = Q)
      
    means, covariances = kf.filter([price_sequence[0]])
    
    filtered_price_sequence = [means[0,0]]
    
    for i in range(1,len(price_sequence)):
        #to track it (streaming)
        new_mean, new_covariance = kf.filter_update(means[-1], covariances[-1], price_sequence[i])
        means = np.vstack([means,new_mean])
        covariances = np.vstack([covariances,new_covariance.reshape((1,3,3))])
    
        filtered_price_sequence.append(means[i,0])
    
    return filtered_price_sequence

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
def main():
    
    #Paramètres du modèle
    #Indiquer à quoi correspond chaque variable
    osigma=0.1;

    transition_matrix = np.array([[1., 0.,0.,0.,0.,0.],[1., 1.,0.,0.,0.,0.],[0.,0,1,0.,0.,0.]])
    transition_covariance = np.zeros((6,6));
    
    

    observation_matrix = np.array([[0., 1.,0.],[0., 0.,1.]])
    observation_covariance = np.eye(2)*osigma

    initial_state_mean = np.array([1,0,10]) # m(0)
    initial_state_covariance = np.eye(3);# p(0)
    
    
    #Observations
    observations = csv.reader(open('voitureObservations.csv','r'),delimiter=' ')
        

    # Filtrage à la main
    # Quels sont les paramètres du constructeur ?
    kf = KalmanFilter(
        transition_matrix, observation_matrix,
        transition_covariance, observation_covariance,
    )    
    
    # Que conserverons les variables suivantes ?
    hand_state_estimates=[initial_state_mean]
    hand_state_cov_estimates=[initial_state_covariance]
    
    for anObs in observations:
        # Quelles étapes du filtrage sont réalisées par la ligne suivante
        (aMean,aCov) = kf.filter_update(hand_state_estimates[-1],hand_state_cov_estimates[-1],anObs)
        hand_state_estimates.append(aMean)
        hand_state_cov_estimates.append(aCov)
    hand_state_estimates=np.array(hand_state_estimates)
    
    # A quoi sert la ligne suivante ?
    hand_positions=np.dot(hand_state_estimates,observation_matrix.T ) # obtenir les observation (H*Mk-1 dans le cours, rouge) projet mon espace etat ver mon espace observation
    
    plt.figure(1)
    plt.plot(observations[:,0],observations[:,1], 'r+')
    plt.plot(hand_positions[:,0],hand_positions[:,1], 'b')
    plt.savefig('illustration_filtrage_main_voiture.pdf')
    plt.close()


    # Filtrage complet
    # Que fait cette section ?
    # Quels sont les paramètres supplémentaires donnés au constructeur ?
    # Quels sont les autres paramètres possibles ?
    kf = KalmanFilter(
        transition_matrix, observation_matrix,
        transition_covariance, observation_covariance,
        initial_state_mean=initial_state_mean, initial_state_covariance=initial_state_covariance,
    )

    
    (filtered_state_estimates,filtered_state_cov_estimates) = kf.filter(observations)
    filtered_positions=np.dot(filtered_state_estimates,observation_matrix.T )
    
    plt.figure(1)
    plt.plot(observations[:,0],observations[:,1], 'r+')
    plt.plot(filtered_positions[:,0],filtered_positions[:,1], 'b')
    plt.savefig('illustration_filtrage_voiture.pdf')
    plt.close()

    # Lissage
    # Que fait cette section ?
    # Quel est la différence avec le filtrage ?
    # Puis-je faire un lissage "à la main" observation par observation ?
   
    kf = KalmanFilter(
        transition_matrix, observation_matrix,
        transition_covariance, observation_covariance,
        initial_state_mean=initial_state_mean, initial_state_covariance=initial_state_covariance,
    )

    (smoothed_state_estimates,smoothed_state_cov_estimates) = kf.smooth(observations) # au lieu de filtre, on appel smooth
    smoothed_positions=np.dot(smoothed_state_estimates,observation_matrix.T )
    
    plt.figure(2)
    plt.plot(observations[:,0],observations[:,1], 'r+')
    plt.plot(smoothed_positions[:,0],smoothed_positions[:,1], 'b')
    plt.savefig('illustration_lissage_voiture.pdf')
    plt.close()

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
n_dim_state = 2#kf.transition_matrix.shape[0]
filtered_state_means = np.zeros((n_timesteps, n_dim_state))
filtered_state_covariances = np.zeros((n_timesteps, n_dim_state, n_dim_state))
#filtered_state_means=measurementsRT
#filtered_state_covariances=measurementsRT
for t in range(n_timesteps - 1):
    measurement=readbgblob()
    measurementsRT[t,:]=measurement
    if t == 0:
        filtered_state_means[t,] = measurements[i,]#kf.initial_state_mean
        filtered_state_covariances[t,] = [[70,70],[70,70]]#kf.initial_state_covariance
    filtered_state_means[t + 1,], filtered_state_covariances[t + 1,] = (
        kf.filter_update(
            filtered_state_means[t,],
            filtered_state_covariances[t,],
            measurementsRT[t,]#,
            #transition_offset=kf.transition_offsets[0], # originally data.transition_offsets[t]
            #observation_offset=kf.observation_offsets
        )
    )
    print filtered_state_means[t+1]
    print measurementsRT[t]

#%% draw estimates
pl.figure()
lines_true = pl.plot(measurementsRT[:,0],measurementsRT[:,1], '-*',color='b')
lines_filt = pl.plot(filtered_state_means[:,0],filtered_state_means[:,1], color='r')
pl.legend((lines_true[0], lines_filt[0]), ('measured', 'filtered'))
pl.show()

#%% sockets for OE

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
# spread_std = np.sqrt(observation_matrix_stepwise.dot(P).dot(observation_matrix_stepwise.transpose())[0][0] + observation_cov)
# print(spread, spread_std)
state_means_stepwise, state_covs_stepwise = kf.filter(observation_stepwise)             # depend on y
# print(state_means_stepwise, state_covs_stepwise)
means_trace.append(state_means_stepwise[0])
covs_trace.append(state_covs_stepwise[0])

for step in range(1, data.shape[0]):
    # print(step)
    x = data[sym_a][step]
    y = data[sym_b][step]
    observation_matrix_stepwise = np.array([[x, 1]])
    observation_stepwise = y

    state_means_stepwise, state_covs_stepwise = kf.filter_update(
        means_trace[-1], covs_trace[-1],
        observation=observation_stepwise,
        observation_matrix=observation_matrix_stepwise)

    # print(state_means_stepwise, state_covs_stepwise)
    # P = covs_trace[-1] + np.eye(2)*state_cov_multiplier                        # This has to be small enough
    # spread = y - observation_matrix_stepwise.dot(means_trace[-1])[0]
    # spread_std = np.sqrt(observation_matrix_stepwise.dot(P).dot(observation_matrix_stepwise.transpose())[0][0] + observation_cov)
    # print(spread, spread_std)
    means_trace.append(state_means_stepwise.data)
    covs_trace.append(state_covs_stepwise)

# ------------------------------------ line evolvement --------------------------------------------#
# colormap
cm = plt.get_cmap('jet')
colors = np.linspace(0.1, 1, data.shape[0])
sc = plt.scatter(data[sym_a], data[sym_b], s=10, c=colors, cmap=cm, edgecolors='k', alpha=0.7)

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
def kalman_track(tser):
    #get the tracking segments
    segments = trigger_and_hold(tser,5,10)
    #delta t
    ts = 0.0002
    #get the first time series slice
    sig_slice = tser[:,0]/np.max(tser[:,0])
    #the kalman filter to smooth in the angle domain
    #array of angle indecies (change to actual dimentions?)
    #set up initial conditions
    #filter for time series
                                                     
    kf2 = KalmanFilter(transition_matrices=np.array([[1.,       ts,       .5*ts**2,     0.,        0,       0, ],
                                                     [0.,       1.,       ts,           0.,        0,       0, ],
                                                     [0.,       0.,       1.,           0.,        0,       0, ],
                                                     [0.,       0.,       0.,           1.,        ts,      0.5*ts**2 ],
                                                     [0.,       0.,       0.,           0.,        1,       ts ],
                                                     [0.,       0.,       0.,           0.,        0 ,      1. ]]),
                  
                  transition_covariance = np.array([[ts**2.,    0.,       0.,           0.,        0,       0.],
                                                     [0.,       1.,       0.,           0.,        0.,      0.],
                                                     [0.,       0,        1/ts**2,      0.,        0.,      0.],
                                                     [0.,       0.,       0.,           ts**2.,    0.,      0.],
                                                     [0.,       0.,       0.,           0.,        1.,      0.],
                                                     [0.,       0,        0.,           0.,        0.,      1/ts**2]]),

                  observation_matrices= np.array([[1,0,0,0,0,0],
                                                  [0,0,0,1,0,0]]),

                  observation_covariance = np.array([[1, 0.,],
                                                     [0.,1]]))
                                                     


    #how much of the time series to filter - for debugging
    testrange = len(segments)-1
    #allocate memory
    filtered_state_means = np.zeros((testrange+1,6))
    filtered_state_cov = np.zeros((testrange+1,6,6))
    
    
    #x = arange(len(sig_slice),dtype = float)
    #kf1 = KalmanFilter(transition_matrices=np.array([[1, 1], [0, 1]]),observation_covariance = [[100000.0]])
    #sm_slice = kf1.smooth(sig_slice)[0][:,0]
    #p = [p_m1, p_m2, sd1, sd2, sc1, sc2]
    #plsq = leastsq(res, p, args = (sm_slice/np.sum(sm_slice), x))
    filtered_state_means[0,:] = [m1,0,0,m2,0,0]
    
    filtered_state_cov[0,:] = np.zeros((6,6))
    #get the valid tracking segments
    
    #itteratively track
    for t in range(testrange):
        print t
        if segments[t]:
            sig_slice = tser[:,t]/np.max(tser[:,t])
            observation = observe(sig_slice,filtered_state_means[t,0],filtered_state_means[t,3])
            filtered_state_means[t+1,:],filtered_state_cov[t+1,:] = kf2.filter_update(filtered_state_means[t,:],
                                                                                      filtered_state_cov[t,:,:],
                                                                                      observation)
        else:
            filtered_state_means[t+1,:] = [m1,0,0,m2,0,0]
            filtered_state_cov[t+1,:] = filtered_state_cov[t,:]
    return filtered_state_means

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
        print ''
        
        #t = 0
        mu = mu[-1]
        sig = sig[-1]
        for i in range(2, len(measurements)):  # starting with third measurement...

            if rospy.is_shutdown():  # easy interruption
                break

            t = i * 1/10.
            #t += 1/10.
            measurement = measurements[i]
            #print measurement
            #print mu
            #print sig
            mu, sig = kf.filter_update(mu, sig, observation=measurement)
            #print mu
            #print sig
            #print ''

            image_temp = image.copy()
            cv2.circle(image_temp, (measurement.data[0], measurement.data[1]), 2, (255,0,0), -1)
            cv2.circle(image_temp, (int(mu[0]), int(mu[1])), 5, (0,0,255))
            spread = numpy.sqrt(sig.diagonal()) * 10
            cv2.rectangle(image_temp, (int(mu[0] - spread[0]), int(mu[1] - spread[1])),
                          (int(mu[0] + spread[0]), int(mu[1] + spread[1])),
                          (0,0,255), 2)
            cv2.imshow('2D Position & Variance', image_temp)
            cv2.waitKey(1000/100)  # display at 100Hz

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
from pykalman import KalmanFilter

kf = KalmanFilter(n_dim_state=1, n_dim_obs=1)
# print y
# output=kf.em(y).smooth(normalizedVals)[0]
##print output
# plot('kalman', output)
# plt.show()


# means, covariances = kf.filter(y[0])
import numpy as np

means = np.zeros((1, 1))
covariances = np.zeros((1, 1))

output = []

next_mean, next_covariance = kf.filter_update(means[-1], covariances[-1], y[0])

for new_measurement in y[1:]:
    next_mean, next_covariance = kf.filter_update(next_mean, next_covariance, new_measurement)

    print next_mean, next_covariance
    output.append(next_mean[0])

print output
plot("Filtered", output)
plt.show()

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
class KalmanFilterPairsTradingStrategy(StrategyBase):
    """
	MeanReversionSpreadStrategy
    """

    def __init__(
            self, events_engine, data_board
    ):
        super(KalmanFilterPairsTradingStrategy, self).__init__(events_engine, data_board)
        self.symbols = ['EWA US Equity', 'EWC US Equity']
        self.bollinger_scaler = 1.0
        self.state_cov_multiplier = np.power(0.01, 2)  # 0.1: spread_std=2.2, cov=16  ==> 0.01: 0.22, 0.16
        self.observation_cov = 0.001

        self.current_ewa_size = 0
        self.current_ewc_size = 0
        self._kf = None                # Kalman Filter
        self._current_state_means = None
        self._current_state_covs = None

    def on_bar(self, bar_event):
        print(bar_event.bar_start_time)
        A = self._data_board.get_hist_price(self.symbols[0], bar_event.bar_start_time)
        B = self._data_board.get_hist_price(self.symbols[1], bar_event.bar_start_time)

        x = A['Price'].iloc[-1]
        y = B['Price'].iloc[-1]
        observation_matrix_stepwise = np.array([[x, 1]])
        observation_stepwise = y
        spread = None
        spread_std = None

        if self._kf is None:
            self._kf = KalmanFilter(n_dim_obs=1, n_dim_state=2,
                              initial_state_mean=np.ones(2),  # initial value
                              initial_state_covariance=np.ones((2, 2)),  # initial value
                              transition_matrices=np.eye(2),  # constant
                              observation_matrices=observation_matrix_stepwise,  # depend on x
                              observation_covariance=self.observation_cov,  # constant
                              transition_covariance=np.eye(2) * self.state_cov_multiplier)  # constant
            P = np.ones((2, 2)) + np.eye(2)*self.state_cov_multiplier
            spread = y - observation_matrix_stepwise.dot(np.ones(2))[0]
            spread_std = np.sqrt(observation_matrix_stepwise.dot(P).dot(observation_matrix_stepwise.transpose())[0][0] + self.observation_cov)
            state_means_stepwise, state_covs_stepwise = self._kf.filter(observation_stepwise)  # depend on y
            self._current_state_means = state_means_stepwise[0]
            self._current_state_covs = state_covs_stepwise[0]
        else:
            state_means_stepwise, state_covs_stepwise = self._kf.filter_update(
                self._current_state_means, self._current_state_covs,
                observation=observation_stepwise,
                observation_matrix=observation_matrix_stepwise)
            P = self._current_state_covs + np.eye(2)*self.state_cov_multiplier                        # This has to be small enough
            spread = y - observation_matrix_stepwise.dot(self._current_state_means)[0]
            spread_std = np.sqrt(observation_matrix_stepwise.dot(P).dot(observation_matrix_stepwise.transpose())[0][0] + self.observation_cov)
            self._current_state_means = state_means_stepwise
            self._current_state_covs = state_covs_stepwise

        # residual is assumed to be N(0, spread_std)
        bollinger_ub = self.bollinger_scaler * spread_std
        bollinger_lb = - self.bollinger_scaler * spread_std
        coeff = self._current_state_means[0]
        print(spread, spread_std)

        if (spread > bollinger_ub) and (self.current_ewa_size >= 0):
            print ('Hit upper band, short spread.')
            o = OrderEvent()
            o.full_symbol = self.symbols[0]
            o.order_type = OrderType.MARKET
            o.order_size = int(-1000*coeff) - self.current_ewa_size
            self.place_order(o)
            self.current_ewa_size = int(-1000*coeff)

            o = OrderEvent()
            o.full_symbol = self.symbols[1]
            o.order_type = OrderType.MARKET
            o.order_size = 1000 - self.current_ewc_size
            self.place_order(o)
            self.current_ewc_size = 1000

        elif (spread < bollinger_lb) and (self.current_ewa_size <= 0):
            print('Hit lower band, long spread.')
            o = OrderEvent()
            o.full_symbol = self.symbols[0]
            o.order_type = OrderType.MARKET
            o.order_size = int(1000 * coeff) - self.current_ewa_size
            self.place_order(o)
            self.current_ewa_size = int(1000 * coeff)

            o = OrderEvent()
            o.full_symbol = self.symbols[1]
            o.order_type = OrderType.MARKET
            o.order_size = -1000 - self.current_ewc_size
            self.place_order(o)
            self.current_ewc_size = -1000

        elif (spread < 0) and (spread > bollinger_lb) and (self.current_ewa_size < 0):
            print('spread crosses below average.cover short spread')
            o = OrderEvent()
            o.full_symbol = self.symbols[0]
            o.order_type = OrderType.MARKET
#.........这里部分代码省略.........

# 需要导入模块: from pykalman import KalmanFilter [as 别名]
# 或者: from pykalman.KalmanFilter import filter_update [as 别名]
def kf_metric_plot(metric):
    assert isinstance(metric, Metric), "Not a metric: {0!s}".format(metric)
    assert hasattr(metric, 'data'), "No Data"

    import numpy as np
    from bounos import DataPackage

    from pykalman import KalmanFilter

    data = DataPackage("/dev/shm/dat-2013-02-01-13-58-48.aietes.npz")

    rnd = np.random.RandomState(0)

    # generate a noisy sine wave to act as our fake observations
    n_timesteps = data.tmax
    x = range(0, n_timesteps)
    records = metric.data

    try:
        obs_dim = len(records[0])
    except TypeError as e:
        obs_dim = 1

    observations = np.ma.array(records)  # to put it as(tmax,3)
    masked = 0
    for i in x:
        try:
            if rnd.normal(2, 2) >= 0:
                observations[i] = np.ma.masked
                masked += 1
        except BaseException as e:
            print(i)
            raise e

    print("{0:f}% Masked".format(((masked * 100.0) / data.tmax)))

    print("Records: Shape: {0!s}, ndim: {1!s}, type: {2!s}".format(records.shape, records.ndim, type(records)))

    # create a Kalman Filter by hinting at the size of the state and observation
    # space.  If you already have good guesses for the initial parameters, put them
    # in here.  The Kalman Filter will try to learn the values of all
    # variables.
    kf = KalmanFilter(n_dim_obs=obs_dim, n_dim_state=obs_dim)

    # You can use the Kalman Filter immediately without fitting, but its estimates
    # may not be as good as if you fit first.

    # states_pred = kf.em(observations, n_iter=data.tmax).smooth(observations)
    # print 'fitted model: %s' % (kf,)

    # Plot lines for the observations without noise, the estimated position of the
    # target before fitting, and the estimated position after fitting.
    fig = plt.figure(figsize=(16, 6))
    ax1 = fig.add_subplot(111)
    filtered_state_means = np.zeros((n_timesteps, kf.n_dim_state))
    filtered_state_covariances = np.zeros(
        (n_timesteps, kf.n_dim_state, kf.n_dim_state))

    for t in x:
        if t == 0:
            tmp = np.zeros(kf.n_dim_state)
            tmp.fill(500)
            filtered_state_means[t] = tmp
            print(filtered_state_means[t])
            continue

        if masked and not observations.mask[t].any():
            ax1.axvline(x=t, linestyle='-', color='r', alpha=0.1)
        try:
            filtered_state_means[t], filtered_state_covariances[t] = \
                kf.filter_update(
                    filtered_state_means[t - 1], filtered_state_covariances[t - 1], observations[t])
        except IndexError as e:
            print(t)
            raise e
    (p, v) = (filtered_state_means, filtered_state_covariances)
    errors = map(np.linalg.norm, p[:] - records[:])

    pred_err = ax1.plot(x, p[:], marker=' ', color='b',
                        label='predictions-x')
    obs_scatter = ax1.plot(x, records, linestyle='-', color='r',
                           label='observations-x', alpha=0.8)
    ax1.set_ylabel(metric.label)
    ax2 = ax1.twinx()
    error_plt = ax2.plot(
        x, errors, linestyle=':', color='g', label="Error Distance")
    ax2.set_yscale('log')
    ax2.set_ylabel("Error")
    lns = pred_err + obs_scatter + error_plt
    labs = [l.get_label() for l in lns]
    plt.legend(lns, labs, loc='upper right')
    plt.xlim(0, 500)
    ax1.set_xlabel('time')

    fig.suptitle(data.title)

    print("Predicted ideal {0!s}: {1!s}".format(metric.label, str(p[-1])))

    plt.show()