% EXAMPLE An example on how to use the Matlab codes CENSOREDML and
% CENSOREDVAR.
%
%   EXAMPLE defines synthetic pathloss data with given pathloss parameters
%   and a given censoring level, c. Then, estimates for the pathloss 
%   parameters of censored data are calculated using ordinary least squares
%   (OLS) and the maximum-likelihood (ML) estimator found in CENSOREDML. 
%   Standard errors for the parameter estimates of the ML estimator are 
%   calculated using CENSOREDVAR. If You use this code, or parts of it, 
%   please cite the following paper: C. Gustafson, T. Abbas, D. Bolin
%   and F. Tufvesson, "Statistical Modeling and Estimation of
%   Censored Pathloss Data", IEEE Wireless Communications Letters, 2015.
%
%   Copyright (C) 2015, David Bolin, Carl Gustafson.
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.

% Define synthetic parameters
d=linspace(10,200,200)'; %Distances [m]
d0=1;                    % Reference distance [m]   
PLd0=20*log10(4*pi*d0/(3e8/5.9e9)); %Reference pathloss value at d0.
n=2;                     %Pathloss exponent.
sigma=4;                 %Large scale standard deviation [dB].
c=90;                    %Censoring level [dB]

% Generate synthetic data 
y=PLd0+10*n*log10(d/d0)+sigma*randn(length(d),1);
x=[ones(length(y),1) 10*log10(d/d0)];

%Define indicator for censore and uncensored samples
t=zeros(1,length(d));
t(find(y<c))=1;
t(find(y>=c))=0;

%Censoring of the synthetic data
yt=zeros(length(d),1);
yt(find(y<c))=y(find(y<c));
yt(find(y>=c))=c;
xt=x(find(y<c),:);

%Ordinary Least Squares (OLS) estimates
a_est = inv(xt'*xt)*xt'*yt(find(t));
sigma2_est = var(yt(find(t))-xt*a_est);

%Maximum likelihood (ML) estimates
thetahat=censoredml(x,yt,c,t,a_est,sigma2_est)

%Asymptotic variance
Avarhat=censoredvar(x,c,thetahat(1),thetahat(2),thetahat(3));
Avar=censoredvar(x,c,PLd0,n,sigma);

%S=['PL(d0)';'n';'sigma'];
disp('   Estimate    OLS     ML    | True value');
disp(sprintf('   PL(d0)      %.2f   %.2f | %.2f', a_est(1), thetahat(1),PLd0))
disp(sprintf('   n           %.3f   %.3f | %.3f', a_est(2), thetahat(2),n))
disp(sprintf('   sigma       %.3f   %.3f | %.3f', sqrt(sigma2_est), thetahat(3),sigma))
disp('   ---------------------------------');
disp('   Asymptotic standard errors (SE) of ML estimates');
disp('                 SE(true value) SE(estimated value)');
disp(sprintf('   SE(PL(d0))    %.2f           %.2f', sqrt(Avar(1)), sqrt(Avarhat(1))))
disp(sprintf('   SE(n)         %.3f          %.3f', sqrt(Avar(2)), sqrt(Avarhat(2))))
disp(sprintf('   SE(sigma^2)   %.3f          %.3f', sqrt(Avar(3)), sqrt(Avarhat(3))))