In [8]:

    

# Generate exact using KF
d = 400
tauPhi = 10.
filename = 'simulatedData/d'+str(d)+'tauPhi'+str(tauPhi)+'y.txt'
y = loadtxt(filename)
filename = 'simulatedData/d'+str(d)+'tauPhi'+str(tauPhi)+'P.txt'
P = loadtxt(filename)
P = P[:d,:d]

T = y.shape[1]
a = 0.5
tauPsi = 1.
tauRho = 1.

Q = P
R = (1/tauPhi)*eye(d)
A = tauRho*a*P

xpred = zeros( (T,d) )
xfilt = zeros( (T,d) )
Ppred = zeros( (T, d, d) )
Pfilt = zeros( (T, d, d) )

for t in range(T):
    #print t
    xpred[t,:] = dot(A,xfilt[t-1,:])
    Ppred[t,:,:] = dot(dot(A,Pfilt[t-1,:,:]),A.T) + Q
    
    S = Ppred[t,:,:] + R
    K = dot(Ppred[t,:,:],inv(S))
    
    xfilt[t,:] = xpred[t,:] + dot(K,y[:,t]-xpred[t,:])
    Pfilt[t,:,:] = Ppred[t,:,:] - dot(K,Ppred[t,:,:])

margP = zeros((T,d))
for t in range(T):
    margP[t,:] = diagonal(Pfilt[t,:,:])
savetxt('exact/d'+str(d)+'tauPhi'+str(tauPhi)+'xfilt.txt',xfilt)
savetxt('exact/d'+str(d)+'tauPhi'+str(tauPhi)+'margPfilt.txt',margP)


    

In [ ]: