Я застрял с этой проблемой. Я пытался реализовать этот алгоритм в python. Но по какой-то причине это не работает. Это в основном следует за наблюдениями, как видно на выходных графиках. Я не могу понять, в чем ошибка модели.
Я использовал описанные здесь дистрибутивы (все гауссовские): https://stats.stackexchange.com/questions/237468/bootstrap-filter-particle-filter-algorithmunderstanding
Алгоритм доступен на страница 9 здесь: https://link.springer.com/content/pdf/10.1007%2F978-1-4757-3437-9_1.pdf
Вывод моего кода доступен здесь, где, как видно, прогнозы просто следуют наблюдениям: https://i.stack.imgur.com/Q3r4R.png
Мой код (для bootstrap повторной выборки последнего состояния достаточно по сравнению со всем путем):
# Reference https://stats.stackexchange.com/questions/237468/bootstrap-filter-particle-filter-algorithmunderstanding
y_sigma = 1
f_sigma = 1
mu_sigma = 1
mu_mean = 0
N= 10000 # number of samples
start_t = 1
end_t = 100
step_t = 1
# generating fake data
n = int(int(end_t-start_t)/step_t + 1) # number os states
X_true = np.cumsum(np.random.normal(0,1,n)) # true (unobservable) values for x
y_obs = np.random.normal(X_true,1,n) # observed values
# initialization, t=0
Sample = np.zeros((N,1))
X_0 = np.random.normal(mu_mean,mu_sigma,N) # N samples from mu ~ Normal(mu_mean,mu_sigma)
X_0 = X_0.reshape(N,1)
Sample[:,0] = np.array([X_0[:,0]]) # an array storing all the samples
X_old = X_0
for t in range(start_t,end_t+step_t,step_t):
f_mean = X_old.reshape(N)
X_new = np.random.normal(f_mean,f_sigma,N)
Sample = np.append(Sample,X_new.reshape(N,1),axis=1)
# importance sampling step
wnew = sp.stats.norm.pdf(y_obs[t-1], X_new, y_sigma) #y_pdf(y_sigma,X_new,y_obs[t-1])
Wnew = wnew / sum(wnew) # normalizing the weights
# selection step
Sample[:,t] = np.random.choice(Sample[:,t], N, p=Wnew)
# updating
X_old = X_new
pred = np.sum(Sample,axis=0)/N
# visualizations
plt.figure(figsize=(20,8))
plt.plot(list(np.arange(start_t,end_t+step_t,step_t)),X_true, color='red', linewidth = 1, label = 'True X')
plt.scatter(list(np.arange(start_t,end_t+step_t,step_t)),y_obs, color='blue', label = 'observations')
plt.plot(list(np.arange(start_t,end_t+step_t,step_t)),pred[1:], color='green', linewidth = 1, label = 'predictions')
print('Error between true X and observations is:', np.sum(abs(y_obs-X_true))/N)
print('Error between true X and predicted values is:', np.sum(abs(pred[1:]-X_true))/N)
plt.legend()
Я также проверил очень простой случай, как показано ниже. Похоже, это просто следствие наблюдений. Например, я могу использовать эти истинные данные и данные наблюдений:
X_true = np.cumsum(np.random.normal(0,1,n)) # true (unobservable) values for x
y_obs = 0*np.random.normal(X_true,1,n)+10 # observed values