%   [x,z]=ExcProbNPU(opt,xd,xu,dx,y,Mmin,lamb,eps,h,xx,ambd)
%
%USING THE NONPARAMETRIC ADAPTATIVE KERNEL APPROACH EVALUATES THE  
%   EXCEEDANCE PROBABILITY VALUES FOR THE UNBOUNDED NONPARAMETRIC 
%   DISTRIBUTION FOR MAGNITUDE. 
%
% AUTHOR: S. Lasocki 06/2014 within IS-EPOS project.
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative 
% to estimating the magnitude distribution functions. It is assumed that 
% the magnitude distribution is unlimited from the right hand side. 
% The estimation makes use of the previously estimated parameters of kernel 
% estimation, namely the smoothing factor, the background sample and the 
% scaling factors for the background sample. The background sample 
% - xx comprises the randomized values of observed magnitude doubled 
% symmetrically with respect to the value Mmin-eps/2.
% The exceedance probability of magnitude M' in the time period of length 
% T' is the probability of an earthquake of magnitude M' or greater to 
% occur in T'. 
% Depending on the value of the parameter opt the exceedance probability 
% values are calculated for a fixed time period T' and different magnitude 
% values or for a fixed magnitude M' and different time period length 
% values. In either case the independent variable vector starts from  
% xd, up to xu with step dx. In either case the result is returned in the
% vector z.
%
% REFERENCES:
%Silverman B.W. (1986) Density Estimation fro Statistics and Data Analysis, 
%   Chapman and Hall, London 
%Kijko A., Lasocki S., Graham G. (2001) Pure appl. geophys. 158, 1655-1665
%Lasocki S., Orlecka-Sikora B. (2008) Tectonophysics 456, 28-37
%       
% INPUT:
%   opt - determines the mode of calculations. opt=0 - fixed time period
%       length (y), different magnitude values (x), opt=1 - fixed magnitude 
%       (y), different time period lengths (x)
%   xd - starting value of the changeable independent variable 
%   xu - ending value of the changeable independent variable
%   dx - step change of the changeable independent variable
%   y - fixed independent variable value: time period length T' if opt=0, 
%        magnitude M' if opt=1
%   Mmin - lower bound of the distribution - catalog completeness level
%   lamb - mean activity rate for events M>=Mmin
%   eps - length of the round-off interval of magnitudes.
%   h - kernel smoothing factor.
%   xx - the background sample
%   ambd - the weigthing factors for the adaptive kernel
%
% OUTPUT:
%   x - vector of changeable independent variable: magnitudes if opt=0, 
%       time period lengths if opt=1, 
%       x=(xd:dx:xu)
%   z - vector of exceedance probability values of the same length as x
%
% LICENSE
%     This file is a part of the IS-EPOS e-PLATFORM.
%
%     This 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.
% 

function [x,z]=ExcProbNPU(opt,xd,xu,dx,y,Mmin,lamb,eps,h,xx,ambd)

% -------------- VALIDATION RULES -------------  K_21NOV2016
 if dx<=0;error('Step must be greater than 0');end
%----------------------------------------------------------


x=(xd:dx:xu)';
n=length(x);

if opt==0
    for i=1:n
        CDF_NPU=2*(Dystr_npr(x(i),xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
        z(i)=1-exp(-lamb*y.*(1-CDF_NPU));
    end
else
    CDF_NPU=2*(Dystr_npr(y,xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
    z=1-exp(-lamb*(1-CDF_NPU).*x);
end
end


function [Fgau]=Dystr_npr(y,x,ambd,h)

%Nonparametric adaptive cumulative distribution for a variable from the
%interval (-inf,inf)

% x - the sample data 
% ambd - the local scaling factors for the adaptive estimation 
% h - the optimal smoothing factor 
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)

n=length(x);
Fgau=sum(normcdf(((y-x)./ambd')./h))/n;
end

