***************LED_data************************
LED_data={{T,groupI},{T,groupII}};
***************Main functions************************
(1) IGWienerNP--The EM algorithm for the semiparametric Wiener process model with IG subordinator.
function out=IGWienerNP(Obs,opt)
% 
% Obs: the degradation observations, an n-length array
%         each element of Obs corresponds to the observation of one group
% Obs{i}: a length-2 array, the first element is an Ki*1 vector 
%              representing the observation times 
%              the second element is an Ki*Mi matrix representing the degradation
%              observations, where each column is the degradation of one unit.
% opt: if opt=='disp', the detailed estimation result will be displayed.
%
% out={T,Lam,sigma2,zeta,loglike}
%
% Example: IGWienerNP(LED_data);
%
(2) IGWienerP--The EM algorithm for the parametric Wiener process model with IG subordinator.
function out=IGWienerP(Obs,Lam,opt)
% Obs: the degradation observations, an n-length array
%         each element of Obs corresponds to the observation of one group
% Obs{i}: a length-2 array, the first element is an Ki*1 vector 
%              representing the observation times 
%              the second element is an Ki*Mi matrix representing the degradation
%              observations, where each column is the degradation of one unit.
% Lam: is a structure corresponding to the definition Lam(t)=beta*S(t):
%         which contains .funcS and .funcdS_da (both are function handles)
% For example, corresponding to Lam(t)=beta*t^alpha, we should define
%     funcS=@(t,a) t.^a;
%     funcdS_da=@(t,a) t.^a.*log(t);
%     Lam=struct('funcS',funcS,'funcdS_da',funcdS_da);
%     If Lam is not defined, by default, we have Lam=beta*t, i.e., the expected degradation path is linear
% opt: if opt=='disp', the detailed estimation result will be displayed.
%
% out={beta,alpha,sigma2,zeta}
%
% Example: 
funcS=@(t,a) t.^a;
funcdS_da=@(t,a) t.^a.*log(t);
Lam=struct('funcS',funcS,'funcdS_da',funcdS_da);
IGWienerP(LED_data, Lam)
%
(3) IGWienerP_infoM--The observed information for the parametric Wiener process model with IG subordinator.
[scoreV,infoM]=IGWienerP_infoM(Obs,estP, Lam)
% Obs: the degradation observations, an n-length array
%         each element of Obs corresponds to the observation of one group
% Obs{i}: a length-2 array, the first element is an Ki*1 vector 
%              representing the observation times 
%              the second element is an Ki*Mi matrix representing the degradation
%              observations, where each column is the degradation of one unit.
% estP: define the point estimates for the model parameters in the parametric model with Lam(t)=beta*t^alpha
%         i.e., estP=struct('beta',estimated_beta, 'alpha',estimated_alpha, 'sigma2',estimated_sigma2, 'zeta',estimated_zeta);
%        OR estP=struct('beta',estimated_beta, 'sigma2',estimated_sigma2, 'zeta',estimated_zeta,);
%        If alpha is not defined, the defaut linear degradation trend Lam(t)=beta*t will be used.
% Lam: is a structure corresponding to the definition Lam(t)=beta*S(t):
%         which contains .funcS and .funcdS_da (both are function handles)
%         For example, corresponding to Lam(t)=beta*t^alpha, we should define
%         funcS=@(t,a) t.^a;
%         funcdS_da=@(t,a) t.^a.*log(t);
%         funcddS_da2=@(t,a) (t.^a).*(log(t).^2);
%         Lam=struct('funcS',funcS,,'funcdS_da',funcdS_da,'funcddS_da2',funcddS_da2);
%
% scoreV and infoM are the score vector and the information matrix at the estimated [beta,sigma2,zeta,alpha(if applicable)] with the observation Obs
% 
% Example: 
funcS=@(t,a) t.^a;
funcdS_da=@(t,a) t.^a.*log(t);
funcddS_da2=@(t,a) (t.^a).*(log(t).^2);
Lam=struct('funcS',funcS,'funcdS_da',funcdS_da,'funcddS_da2',funcddS_da2);
res=IGWienerP(LED_data, Lam);
estP=struct('beta',res(1), 'alpha',res(2), 'sigma2',res(3), 'zeta',res(4));
[a,b]=IGWienerP_infoM(LED_data,estP,Lam)
%
%
(4) IGWienerP_Tp--The point estimate and asymptotic variance for the p-percent lifetime under the parametric Wiener process model with IG subordinator.
[Tp,var_Tp]=IGWienerP_Tp(estP,estCov,d,p, Lam)
% The quantile T and its variance at a given p based on delta method
%  estP: define the point estimates for the model parameters in the parametric model with Lam(t)=beta*t^alpha
%           i.e., estP=struct('beta',estimated_beta, 'sigma2',estimated_sigma2, 'zeta',estimated_zeta, 'alpha',estimated_alpha);
%           OR estP=struct('beta',estimated_beta, 'sigma2',estimated_sigma2, 'zeta',estimated_zeta,);
%           If alpha is not defined, the defaut linear degradation trend Lam(t)=beta*t will be used.
% estCov: the asymptotic variance at estP
% d: failure threshold
% p: define the p-percent lifetime
% Lam=struct('funcS',funcS,'funcdS_da',funcdS_da,'funcdS_dt',funcdS_dt);
%
% Tp and var_Tp are the point estimate and the asymptotic variance for the p-percent lifetime
%
% Example: 
funcS=@(t,a) t.^a;
funcdS_da=@(t,a) t.^a.*log(t);
funcddS_da2=@(t,a) (t.^a).*(log(t).^2);
Lam=struct('funcS',funcS,'funcdS_da',funcdS_da,'funcddS_da2',funcddS_da2);
res=IGWienerP(LED_data, Lam);
estP=struct('beta',res(1), 'alpha',res(2), 'sigma2',res(3), 'zeta',res(4));
[a,b]=IGWienerP_infoM(LED_data,estP,Lam);
Lam.funcdS_dt=@(t,a) a*t.^(a-1);
[Tp,var_Tp]=IGWienerP_Tp(estP,inv(b),0.15,0.5,Lam);
%

***************Other functions************************
(6) IGWienerNP_lnL: give the loglikelihood value for the parametric Wiener process model with IG subordinator.
loglike=IGWienerNP_lnL(Obs,estP)
% Obs: the degradation observations, an n-length array as before.
% estP: define the point estimates for the model parameters in the semiparametric model
% i.e., estP=estpara=struct('zeta',estimated_zeta,'sigma2',estimated_sigma2,'dLam',estimated_delta_Lam_at_the_superset_of_T);
%
%
(5) IGWienerP_lnL: give the loglikelihood value for the parametric Wiener process model with IG subordinator.
loglike=IGWienerP_lnL(Obs,estP,funcS)
% Obs: the degradation observations, an n-length array as before.
% estP: define the point estimates for the model parameters in the parametric model with Lam(t)=beta*t^alpha
% i.e., estP=struct('beta',estimated_beta, 'sigma2',estimated_sigma2, 'zeta',estimated_zeta, 'alpha',estimated_alpha);
% OR estP=struct('beta',estimated_beta, 'sigma2',estimated_sigma2, 'zeta',estimated_zeta,);
%
% funcS: the function handle for S(t;alpha)
