My Particle filter is not providing approx parameters of the model

Question

Hellow all,

Particle Filter

I am trying to estimate parameters of a model using particle filter. I am using simple particle filter to estimate parameters of model. But estimate parameters are not approx near to true parameters and iterations are not converging.

Can anyone help me optimize the particle filter code to get better estimate parameters and converge the iteration?

I have provided code of particle filter for your reference.

# Universal Failure Rate Function
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Y = 0         # Base
K = 10       # Scale
beta = 4    # Shape
alpha = 14    # Scale
num = 500
measurement_noise_std = .3

time = np.linspace(0, 10, num)
UFRF = Y + K*(beta/alpha**beta)*time**(beta-1)

#Add noise to the true data
MFRF = UFRF + np.random.normal(0, measurement_noise_std, size=num)
#MFRF = UFRF + 0.7 * np.random.normal(size=len(UFRF))
plt.plot(time, MFRF)
plt.plot(time, UFRF)

np.random.seed(42) 
beta_space = (2, 8)
alpha_space = (8, 15)
Y_base_space = (-2, 2)
K_scale_space = (5, 15)
num_particles = 1000
motion_model_std = 0.7
std_dev = 0.4

Step 1: Initialization

def initialize_particles(num_particles, Y_base_space, K_scale_space, beta_space, alpha_space):
   Y_bases = np.random.uniform(Y_base_space[0], Y_base_space[1], size=num_particles)
   K_scales = np.random.uniform(K_scale_space[0], K_scale_space[1], size=num_particles)
   betas = np.random.uniform(beta_space[0], beta_space[1], size=num_particles)
   alphas = np.random.uniform(alpha_space[0], alpha_space[1], size=num_particles)
   weights = np.ones(num_particles) / num_particles
   return Y_bases, K_scales, betas, alphas, weights

Step 2: Prediction

def predict_particles(Y_bases, K_scales, betas, alphas, motion_model_std):
   Y_bases = Y_bases + np.random.normal(0, motion_model_std, len(Y_bases))
   K_scales = K_scales + np.random.normal(0, motion_model_std, len(K_scales))
   betas = betas + np.random.normal(0, motion_model_std, len(betas))
   alphas = alphas + np.random.normal(0, motion_model_std, len(alphas))
return Y_bases, K_scales, betas, alphas

Step 3: Update

def update_weight(observed_data, predicted_data, motion_model_std):
   residual = observed_data - predicted_data
   likelihood = np.exp(-0.5 * np.sum(residual**2) / motion_model_std**2)
   #print(likelihood)
   # Handle potential division by zero
   if np.isclose(likelihood, 0):
      return 1e-100  # Small epsilon value
   else:
return likelihood   

def neff(weights):
return 1. / np.sum(np.square(weights))

Step 4: Resampling

def resample(Y_bases, K_scales, betas, alphas, weights):
   indices = np.random.choice(len(Y_bases), size=len(Y_bases), p=weights)
   Y_bases = Y_bases[indices]
   K_scales = K_scales[indices]   
   betas = betas[indices]
   alphas = alphas[indices] 
   N = len(indices)
   if neff(weights) < N:
      #print("true")
      weights = np.ones(len(weights)) / len(weights)  # Reset weights
   else:
     #print("false")
     weights = weights
return Y_bases, K_scales, betas, alphas, weights    # Reset weights

Main Particle Filter Algorithm for Universal Universal Failure Rate Function Parameter Estimation

def particle_filter_MFRF_parameter_estimation(num_particles, Y_base_space, K_scale_space, beta_space, alpha_space, motion_model_std, MFRF, time):

  MFRF = pd.DataFrame(MFRF)  # dataframe
  time = pd.DataFrame(time)  # dataframe

  # Initialization
  Y_bases, K_scales, betas, alphas, weights = initialize_particles(num_particles, Y_base_space, K_scale_space, 
                                                                  beta_space, alpha_space)

  # Estimate
  Y_base_estimate = np.sum(Y_bases * weights)
  K_scale_estimate = np.sum(K_scales * weights)
  beta_estimate = np.sum(betas * weights)
  alpha_estimate = np.sum(alphas * weights)

  predicted_curve = Y_base_estimate + K_scale_estimate * (beta_estimate/alpha_estimate**
                                                            beta_estimate)*time**(beta_estimate-1)
  residual = MFRF - predicted_curve
  std = pd.Series(np.std(residual))
  std = std.values[0]
  print(std)
  while std >= std_dev:
    
      # Prediction
      Y_bases, K_scales, betas, alphas = predict_particles(Y_bases, K_scales, betas, alphas, motion_model_std)
       
      # Update
      for i in range(num_particles):
        
          # Y + K*(beta/eta**beta)*time**(beta-1)
          predicted_MFRF = Y_bases[i] + K_scales[i] * (betas[i]/alphas[i]**betas[i])*time**(betas[i]-1)
          weights[i] = update_weight(MFRF, predicted_MFRF, motion_model_std)
        
      weights /= np.sum(weights)

      # Resampling
      Y_bases, K_scales, betas, alphas, weights = resample(Y_bases, K_scales, betas, alphas, weights)
    
    Y_base_estimate = np.sum(Y_bases * weights)
    K_scale_estimate = np.sum(K_scales * weights)
    beta_estimate = np.sum(betas * weights)
    alpha_estimate = np.sum(alphas * weights)

    predicted_curve = Y_base_estimate + K_scale_estimate * (beta_estimate/alpha_estimate**
                                                            beta_estimate)*time**(beta_estimate-1)
    residual = MFRF - predicted_curve
    std = pd.Series(np.std(residual))
    std = std.values[0]
    print(std)
 
# Estimate the parameters as the weighted mean of particles
plt.plot(Y_bases)
plt.plot(K_scales)
plt.plot(betas)
plt.plot(alphas)

return Y_base_estimate, K_scale_estimate, beta_estimate, alpha_estimate

Main Function Execution

Y_base_estimate, K_scale_estimate, beta_estimate, alpha_estimate = particle_filter_MFRF_parameter_estimation(
num_particles, Y_base_space, K_scale_space, beta_space, alpha_space, motion_model_std, MFRF, time)

Plot

predicted_MFRF = Y_base_estimate + K_scale_estimate * (beta_estimate/alpha_estimate**beta_estimate)*time**(beta_estimate-1)

print("Y_base_estimate = ",Y_base_estimate)
print("K_scale_estimate = ", K_scale_estimate)
print("beta_estimate = ", beta_estimate)
print("alpha_estimate = ", alpha_estimate)
plt.plot(time, predicted_MFRF)
plt.plot(time, MFRF)