function [thetaHat,rul]=BM_Battery(para0,weigh)
global DegraUnit initDisPar ...
TimeUnit time y thres ParamName thetaTrue signiLevel ns ny np
%=== PROBLEM DEFINITION 1 (Required Variables) ==============
WorkName='Battery_BM';  % work results are saved by WorkName
DegraUnit='C/1 Capacity';                  % degradation unit
TimeUnit='Cycles';          % time unit (Cycles, Weeks, etc.)
time=[0:5:200]';    % time at both measurement and prediction
y=[1.00 0.99 0.99 0.94 0.95 0.94 0.91 0.91 0.87 0.86]'; %[nyx1]: measured data
thres=0.7;                       % threshold (critical value)
ParamName=['b'; 's'];%[npx1]: parameters' name to be estimated
initDisPar=[0 0.02; 1e-5 0.1]; %[npx2]: prob. parameters of init./prior dist
thetaTrue=[0.003; 0.02];   %[npx1]: true values of parameters
signiLevel=5;          % significance level for C.I. and P.I.
ns=5e3;                         % number of particles/samples
burnIn=0.2;                                 % ratio for burn-in
%============================================================
% % % PROGNOSIS using BM with MCMC
ny=length(y);
np=size(ParamName,1);
sampl(:,1)=para0;                          %% Initial Samples
[~,jPdf0]=MODEL(para0,time(1:ny));     %% Initial (joint) PDF
for i=2:ns/(1-burnIn);                        %% MCMC Process
 % proposal distribution (uniform)
 para1(:,1)=unifrnd(para0-weigh,para0+weigh);
 % (joint) PDF at new samples
 [~,jPdf1]=MODEL(para1,time(1:ny));
 % acception/rejection criterion
 if rand<min(1,jPdf1/jPdf0)
  para0=para1; jPdf0=jPdf1;
 end;
 sampl(:,i)=para0;
end;
nBurn=ns/(1-burnIn)-ns;                     %% Sampling Trace
figure(4);for j=1:np;subplot(np,1,j);plot(sampl(j,:));hold on
plot([nBurn nBurn],[min(sampl(j,:)) max(sampl(j,:))],'r');end
thetaHat=sampl(:,nBurn+1:end);      %% Final Sampling Results
for k=1:length(time(ny:end));       %% Degradation Prediction
 [zHat(k,:),~]=MODEL(thetaHat,time(ny-1+k));
 degraPredi(k,:)=normrnd(zHat(k,:),thetaHat(end,:));
end;
% % % POST-PROCESSING
degraTrue=[];                             %% True Degradation
if ~isempty(thetaTrue); degraTrue=MODEL(thetaTrue,time); end
rul=POST(thetaHat,degraPredi,degraTrue); %% RUL & Result Disp
Name=[WorkName ' at ' num2str(time(ny)) '.mat']; save(Name);
end
function [z,poste]=MODEL(param,t)
 global y ParamName initDisPar ny np
 for j=1:np; eval([ParamName(j,:) '=param(j,:);']); end
 %===== PROBLEM DEFINITION 2 (model equation) ===============
 z=exp(-b.*t);
 %===========================================================
 if length(y)~=length(t); poste=[];
 else                                                 
  prior= ...                                         %% Prior
        prod(unifpdf(param,initDisPar(:,1),initDisPar(:,2)));
  likel=(1./(sqrt(2.*pi).*s)).^ny.* ...         %% Likelihood
          exp(-0.5./s.^2.*(y-z)'*(y-z));
  poste=likel.*prior;                            %% Posterior
 end
end