Particle filter error solution

Hi, I am still new to malab. Try to run a simple particle filter but error keep appear saying too many input arguments and error in

(line 77 ) [xh(:,k), pf] = particle_filter(xk, Bond, y(:,k), pf,'systematic_resampling'); saying unrecognized y function. Please help me.

%% clear memory, screen, and close all figures
clear, clc, close all;
%% initial setting
T = 4;     %step number
Ns = 10;   %number of ensemble
MB = 60;    %number of bond to 100% cured
%% Predictive ensemble formation
nx  = 3;    %numberofstate (k, Bond, vk)
%nx  = size(xk1);
xk = zeros(233,5800,67);   %taktahubetulke tak
for j = 1:Ns
end
Bond  = xk(3,:,:);
Fbond = zeros(Ns,1);
for j = 1:Ns
    Fbond(j,:) = Bond(j);
end
%% Observed value matrix
ny = 1;    %number of observations
Y = '(0 ; 16.7)';
%yk  = @(k) Y(ny,k);
yk  = Y(ny,2);
%yek = zeros(ny,Ns,m);
%yek = zeros(ny,Ns);
%yek(ny,1:Ns) = yk; 
%obs  = @(k, xk, vk) (1-(xk(3,1:Ns)/MB)) + vk;
obs  = @(k, Bond, vk) ((Bond/MB)*100) + vk;
%% process noise and noise generator function
nu = 1;                                           % size of the vector of process noise
sigma_u = sqrt(10);
p_sys_noise   = @(u) normpdf(u, 0, sigma_u);
gen_sys_noise = @(u) normrnd(0, sigma_u);         % sample from p_sys_noise (returns column vector)
%% observation noise and noise generator function
nv = 1;                                           % size of the vector of observation noise
sigma_v = sqrt(1);
p_obs_noise   = @(v) normpdf(v, 0, sigma_v);
gen_obs_noise = @(v) normrnd(0, sigma_v);         % sample from p_obs_noise (returns column vector)
%% Likelihood p(y[k] | x[k])
% (under the suposition of additive process noise)
%p_yk_given_xk = @(k, yk, xk) p_obs_noise(yk - obs(k, xk, 0));
R = cov(Fbond);
test = (yk-obs(2,Bond,0));
p_yk_given_xk = @(k, yk, Bond) (1/((2*pi)*det(R))^(1/2))*exp((-1/2)*(yk-obs(k,Bond,0))*(R^(-1))*(yk-obs(k,Bond,0)));
%p_yk_given_xk = @(k, yk, Bond) (1/((2*pi)*det(R))^(1/2))*exp((-1/2)*(test)'*(R^(-1))*(test));
YX = (1/((2*pi)*det(R))^(1/2))*exp((-1/2)*(yk-obs(2,Bond,0)).*(R^(-1)).*(yk-obs(2,Bond,0)));
%% %% Separate memory space
x = zeros(nx,T);  y = zeros(ny,T);
u = zeros(nu,T);  v = zeros(nv,T);
%% Separate memory
xh0 = 0;
xh = zeros(nx, T); xh(:,1) = xh0;
yh = zeros(Ns, T); yh(:,1) = obs(1, xh0, 0);
gen_x0 = @(x) normrnd(0, sqrt(10));
pf.k               = 1;                   % initial iteration number
pf.Ns              = 10;                 % number of particles
pf.w               = zeros(pf.Ns, T);     % weights
pf.particles       = zeros(nx, pf.Ns, T); % particles
pf.gen_x0          = gen_x0;              % function for sampling from initial pdf p_x0
pf.p_yk_given_xk   = p_yk_given_xk;       % function of the observation likelihood PDF p(y[k] | x[k])
pf.gen_sys_noise   = gen_sys_noise;       % function for generating system noise
%pf.p_x0 = p_x0;                          % initial prior PDF p(x[0])
%pf.p_xk_given_ xkm1 = p_xk_given_xkm1;   % transition prior PDF p(x[k] | x[k-1])
%% Estimate state
for k = 2:T
    fprintf('Iteration = %d/%d\n',k,T);
    % state estimation
    pf.k = k;
    %[xh(:,k), pf] = particle_filter(sys, y(:,k), pf, 'multinomial_resampling');
    [xh(:,k), pf] = particle_filter(xk, Bond, y(:,k), pf,'systematic_resampling');   
    
    % filtered observation
    yh(:,k) = obs(k, Bond(:,k), 0);
    
end
%% plot of the state vs estimated state by the particle filter vs particle paths
figure
hold on;
h1 = plot(1:T,squeeze(pf.particles),'y');
h2 = plot(1:T,x,'b','LineWidth',4);
h3 = plot(1:T,xh,'r','LineWidth',4);
h4 = plot(1:T,xhmode,'g','LineWidth',4);
legend([h2 h3 h4 h1(1)],'state','mean of estimated state','mode of estimated state','particle paths');
title('State vs estimated state by the particle filter vs particle paths','FontSize',14);

Answer 1

[Walter Roberson]

My guess is that you are using the particle_filter function from https://www.mathworks.com/matlabcentral/fileexchange/35468-particle-filter-tutorial . That function expects sys, yk, pf, resampling_strategy where sys is an aonymous function, pf is a struct, and that is the correct spelling of the resampling_strategy option.

The function does not support passing in a pair of numeric arrays.

Look at the commented out line directly before it to see what it expects.

It is possible, of course, that you are expecting to use a modified version of that code... if so then you need to have that modified version on your path (before the original version if the original version is still present.)

[シティヌルシュハダ モハマド ナシル]

Thank you for your answer, I really appreciate it. As I am still new to matlab, there some word that I don't understand. Do you have other recommendation for me to do with the coding and which part that I need to edit to solve the error.

And how to make sure the modified version is on my path if the original version still present?

[Walter Roberson]

xk = zeros(233,5800,67);   %taktahubetulke tak

You initialize xk to all zero, which is possibly a good thing to do.

for j = 1:Ns
end

You have a loop there, but you do not do anything inside the loop. It is not clear why you have that loop.

Bond  = xk(3,:,:);

xk has not changed since it was initialized, so xk is still all zero. It is not clear why you would pick out the third row specifically out of the all-zero array.

Fbond = zeros(Ns,1);
for j = 1:Ns
    Fbond(j,:) = Bond(j);
end

Bond is all zero because it was taken from the all-zero xk, so your loop is overwriting the Fbond = zeros with values that are all zero. It is not clear why you would do that.

for k = 2:T
    fprintf('Iteration = %d/%d\n',k,T);
    % state estimation
    pf.k = k;
    %[xh(:,k), pf] = particle_filter(sys, y(:,k), pf, 'multinomial_resampling');
    [xh(:,k), pf] = particle_filter(xk, Bond, y(:,k), pf,'systematic_resampling');   
    
    % filtered observation
    yh(:,k) = obs(k, Bond(:,k), 0);
    
end

You have not changed xk since you assigned it all zero, and you do not change it inside the loop. You have not changed Bond since you assigned it all zero, and you do not change it inside the loop. It is not clear why you would want to pass xk or Bond in to particle_filter() instead of the sys that it was designed to work with ?

The particle_filter is designed for the situation in which particle locations evolve according to a function handle, with the function handle named sys for this code. If you want the particle locations to stay constant 0, then instead of editing the particle filter code to try to force it to work with constant arrays instead of with a function handle, then I would suggest you should instead just pass in

sys = @(k, xkm1, uk) zeros(size(xkm1,1),1)

or something like that.

I do not know why you changed the interface to particle_filter so it is difficult to make recommendations.

[シティヌルシュハダ モハマド ナシル]

Thank you very much for your explanation about the coding and suggestion. I really appreciate it.

Actually I was just trying my senior's code without knowing the reason and explanation about the code. My senior is no longer in university and I cannot contact him, so I am very clueless.

[シティヌルシュハダ モハマド ナシル]

Actually my research is about the estimation of curing reaction of epoxy resin using particle filter. I had done simulation and also the experiment about curing reaction and plan to use particle filter to get the accurate curing reaction degree. For the ensemble member, I had done about 50 type of curing reaction simulation with different reaction distance and I understand that this data will act as particle. And the data from the experiment will act as the observation data? (correct me if I am wrong)

I plan to use and modify the particle_filter function https://www.mathworks.com/matlabcentral/fileexchange/35468-particle-filter-tutorial

but I still confuse where I can input all the curing reaction degree data into the coding. I really really appreciate if you can help me. There is no one in my lab that understand the particle filter coding.

%% clear memory, screen, and close all figures
clear, clc, close all;
%% Process equation x[k] = sys(k, x[k-1], u[k]);
nx = 1;  % number of states
sys = @(k, xkm1, uk) xkm1/2 + 25*xkm1/(1+xkm1^2) + 8*cos(1.2*k) + uk;
%% Observation equation y[k] = obs(k, x[k], v[k]);
ny = 1;                                           % number of observations
obs = @(k, xk, vk) xk(1)+ vk;                  % (returns column vector)
%% PDF of process noise and noise generator function
nu = 1;                                           % size of the vector of process noise
sigma_u = sqrt(10);
p_sys_noise   = @(u) normpdf(u, 0, sigma_u);
gen_sys_noise = @(u) normrnd(0, sigma_u);         % sample from p_sys_noise (returns column vector)
%% PDF of observation noise and noise generator function
nv = 1;                                           % size of the vector of observation noise
sigma_v = sqrt(1);
p_obs_noise   = @(v) normpdf(v, 0, sigma_v);
gen_obs_noise = @(v) normrnd(0, sigma_v);         % sample from p_obs_noise (returns column vector)
%% Initial PDF
% p_x0 = @(x) normpdf(x, 0,sqrt(10));             % initial pdf
gen_x0 = @(x) normrnd(0, sqrt(10));               % sample from p_x0 (returns column vector)
%% Transition prior PDF p(x[k] | x[k-1])
% (under the suposition of additive process noise)
% p_xk_given_xkm1 = @(k, xk, xkm1) p_sys_noise(xk - sys(k, xkm1, 0));
%% Observation likelihood PDF p(y[k] | x[k])
% (under the suposition of additive process noise)
p_yk_given_xk = @(k, yk, xk) p_obs_noise(yk - obs(k, xk, 0));
%% Number of time steps
T = 40;
%% Separate memory space
x = zeros(nx,T);  y = zeros(ny,T);
u = zeros(nu,T);  v = zeros(nv,T);
%% Simulate system
xh0 = 0;                                  % initial state (h, Vf)
u(:,1) = 0;                               % initial process noise
v(:,1) = gen_obs_noise(sigma_v);          % initial observation noise
x(:,1) = xh0;
y(:,1) = obs(1, xh0, v(:,1));
for k = 2:T
    % here we are basically sampling from p_xk_given_xkm1 and from p_yk_given_xk
    u(:,k) = gen_sys_noise();              % simulate process noise
    v(:,k) = gen_obs_noise();              % simulate observation noise
    x(:,k) = sys(k, x(:,k-1), u(:,k));   % simulate state
    y(:,k) = obs(k, x(:,k),   v(:,k));     % simulate observation
end
fprintf('Finish simulate system \n')
%% Separate memory
xh = zeros(nx, T); xh(:,1) = xh0;
yh = zeros(ny, T); yh(:,1) = obs(1, xh0, 0);
pf.k               = 1;                   % initial iteration number
pf.Ns              = 10;                 % number of particles
pf.w               = zeros(pf.Ns, T);     % weights
pf.particles       = zeros(nx, pf.Ns, T); % particles
pf.gen_x0          = gen_x0;              % function for sampling from initial pdf p_x0
pf.p_yk_given_xk   = p_yk_given_xk;       % function of the observation likelihood PDF p(y[k] | x[k])
pf.gen_sys_noise   = gen_sys_noise;       % function for generating system noise
%pf.p_x0 = p_x0;                          % initial prior PDF p(x[0])
%pf.p_xk_given_ xkm1 = p_xk_given_xkm1;   % transition prior PDF p(x[k] | x[k-1])
%% Estimate state
for k = 2:T
    fprintf('Iteration = %d/%d\n',k,T);
    % state estimation
    pf.k = k;
    %[xh(:,k), pf] = particle_filter(sys, y(:,k), pf, 'multinomial_resampling');
    [xh(:,k), pf] = particle_filter(sys, y(:,k), pf, 'systematic_resampling');   
    
    % filtered observation
    yh(:,k) = obs(k, xh(:,k), 0);
end
%% Make plots of the evolution of the density
Disp=3;%表示させたい状態変数（1:厚さ，2?6:Vf．6は表面）
figure
hold on;
xi = 1:T;
yi = -25:0.25:25; %変数としてDispを想定した領域に設定
[xx,yy] = meshgrid(xi,yi);
den = zeros(size(xx));
xhmode = zeros(size(xh));
for i = xi
    % for each time step perform a kernel density estimation
    den(:,i) = ksdensity(pf.particles(1,:,i), yi,'kernel','epanechnikov');
    [~, idx] = max(den(:,i));
    % estimate the mode of the density
    xhmode(i) = yi(idx);
    plot3(repmat(xi(i),length(yi),1), yi', den(:,i));
end
view(3);
box on;
title('Evolution of the state density','FontSize',14)
figure
mesh(xx,yy,den);   
title('Evolution of the state density','FontSize',14)
%% plot of the state vs estimated state by the particle filter vs particle paths
figure
hold on;
h1 = plot(1:T,squeeze(pf.particles),'y');
h2 = plot(1:T,x,'b','LineWidth',4);
h3 = plot(1:T,xh,'r','LineWidth',4);
h4 = plot(1:T,xhmode,'g','LineWidth',4);
legend([h2 h3 h4 h1(1)],'state','mean of estimated state','particle paths');
title('State vs estimated state by the particle filter vs particle paths','FontSize',14);
%% plot of the observation vs filtered observation by the particle filter
figure
plot(1:T,y,'b', 1:T,yh,'r');
legend('observation','filtered observation');
title('Observation vs filtered observation by the particle filter','FontSize',14);
return;