%% [NLS]: MATLAB Code for Nonlinear Least Square
 
 
function [thetaHat,rul]=NLS(theta0)
clear global; global DegraUnit 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
y=[ ]';                                %[nyx1]: measured data
thres= ;                         % threshold (critical value)
ParamName=[ ];      %[npx1]: parameters' name to be estimated
thetaTrue=[ ];             %[npx1]: true values of parameters
signiLevel= ;          % significance level for C.I. and P.I.
ns= ;                           % number of particles/samples
%============================================================
% % % PROGNOSIS using NLS
ny=length(y); nt=ny;
np=size(ParamName,1);             %% Deterministic Estimation
Optio=optimset('algorithm',{'levenberg-marquardt',0.01});
[thetaH,~,resid,~,~,~,J]=lsqnonlin(@FUNC,theta0,[],[],Optio);
sse=resid'*resid;                    %% Distribution of Theta
dof=ny-np+1;
sigmaH=sqrt(sse/dof);                 % Estimated std of data
W=eye(ny)*1/sigmaH^2; thetaCov=inv(J'*W*J);
sigTheta=sqrt(diag(thetaCov));       % Estimated std of Theta
mvt=mvtrnd(thetaCov,dof,ns)';               % Generate t-dist
thetaHat=repmat(thetaH,1,ns)+mvt.*repmat(sigTheta,1,ns);
thetaHat(np,:)=sigmaH*ones(1,ns);   %% Final Sampling Results
for k=1:length(time(ny:end));       %% Degradation Prediction
 zHat(k,:)=MODEL(thetaHat,time(ny-1+k));
 degraPredi(k,:)=zHat(k,:)+trnd(dof,1,ns)*sigmaH;
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 objec=FUNC(theta)
 global time y ny
 z=MODEL(theta,time(1:ny));
 objec=(y-z);
end
function z=MODEL(theta,t)
 global ParamName np
 for j=1:np-1; eval([ParamName(j,:) '=theta(j,:);']); end
 %===== PROBLEM DEFINITION 2 (model equation) ===============

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