%% [PF]: MATLAB Code for Particle Filter
 
 
function [thetaHat,rul]=PF
clear global; global DegraUnit initDisPar TimeUnit ...
time y thres ParamName thetaTrue signiLevel ns ny nt np
%=== PROBLEM DEFINITION 1 (Required Variables) ==============
WorkName=' ';            % work results are saved by WorkName
DegraUnit=' ';                             % degradation unit
TimeUnit=' ';               % time unit (Cycles, Weeks, etc.)
time=[ ]';          % time at both measurement and prediction
dt= ;             % time interval for degradation propagation
y=[ ]';                                %[nyx1]: measured data
thres= ;                         % threshold (critical value)
ParamName=[ ];      %[npx1]: parameters' name to be estimated
initDisPar=[ ]; %[npx2]: prob. parameters of init./prior dist
thetaTrue=[ ];             %[npx1]: true values of parameters
signiLevel= ;          % significance level for C.I. and P.I.
ns= ;                           % number of particles/samples
%============================================================
% % % PROGNOSIS using PF
ny=length(y); nt=ny;
np=size(ParamName,1);
for j=1:np;                           %% Initial Distribution
 param(j,:,1)=unifrnd(initDisPar(j,1),initDisPar(j,2),1,ns);
end;
for k=2:length(time);          %% Update Process or Prognosis
 % step1. prediction (prior)
 paramPredi=param(:,:,k-1);
 nk=(time(k)-time(k-1))/dt;
 for k0=1:nk; paramPredi(np,:)=MODEL(paramPredi,dt); end
if k<=ny                                  % (Update Process)
  % step2. update (likelihood)
  likel=normpdf(y(k),paramPredi(np,:),paramPredi(np-1,:));
  % step3. resampling
  cdf=cumsum(likel)./sum(likel);
  for i=1:ns;
   u=rand;
   loca=find(cdf>=u,1); param(:,i,k)=paramPredi(:,loca);
  end;
 else                                           % (Prognosis)
  param(:,:,k)=paramPredi;
 end
end
thetaHat=param(1:np-1,:,ny);        %% Final Sampling Results
paramRearr=permute(param,[3 2 1]);  %% Degradation Prediction
zHat=paramRearr(:,:,np);
degraPredi=normrnd(zHat(ny:end,:),paramRearr(ny:end,:,np-1));
% % % POST-PROCESSING
degraTrue=[];                             %% True Degradation
if ~isempty(thetaTrue); k=1;
 degraTrue0(1)=thetaTrue(np); degraTrue(1)=thetaTrue(np);
 for k0=2:max(time)/dt+1;
  degraTrue0(k0,1)= ...
             MODEL([thetaTrue(1:np-1); degraTrue0(k0-1)],dt);
  loca=find((k0-1)*dt==time,1);
  if ~isempty(loca); k=k+1; degraTrue(k)=degraTrue0(k0);end
 end
end
rul=POST(thetaHat,degraPredi,degraTrue); %% RUL & Result Disp
Name=[WorkName ' at ' num2str(time(ny)) '.mat']; save(Name);
end
function z1=MODEL(param,dt)
 global ParamName np
 for j=1:np; eval([ParamName(j,:) '=param(j,:);']); end
 %===== PROBLEM DEFINITION 2 (model equation) ===============

 %===========================================================
end