%
% CpGFT Matlab file, v.1.0.2011
%
function [status]=CpGFT(type,paths,timepoints,mean_sd,alpha)

%  CpGFT tests if a sample of observed paths of a counting
%  process M(t) fits a given Cox Process(CP), N(t). The hypothesis test assesses
%  the coherence of new observed data with a proposed CP, with known or
%  estimated moments.The authors suggest as an accurate estimation of the
%  intensity the one proposed in [3].
%
%  A goodness-of-fit test for functional data is derived, based on a simple
%  discrete estimator of the intensity process of the CP at each observation 
%  point. The null hypothesis states that the intensity mean is the proposed 
%  one at each observation point. Afterwards, a multiple testing criterion 
%  (Benjamini Hochberg, [4]) yields a decision on the coherence of the observed 
%  complete sample path with the proposed CP (see Bouzas et al. [2]).
%
%  Having observed s sample paths of M(t) of dimension p*m+1, the discrete 
%  estimator of the intensity process l(t) for each sample path is calculated 
%  by splitting it up in m subpaths of p+1 points.The subpath i contains the 
%  values of M(t) at time points t_((i-1)*p) to t_(i*p) for i:1,...,m. Then, 
%  the point estimator proposed in Bouzas et al. [1] is applied. With the s
%  estimated sample paths of the intensity at each time point, an
%  experimental value of the mean of the intensity is calculated and the
%  conclusion of the hypothesis test is achieved at each t_j. Taking into
%  account the information of all isolated tests, the simultaneous
%  inference leads us to a final decision.
%
%  Syntax: function [status]=CpGFT(type,paths,timepoints,mean_sd,alpha)
%
%  Inputs:
%      type         - 0 if the intensity process is unknown and its
%                     properties has been estimated by the user in an 
%                     external procedure to obtain its mean and its
%                     standard deviation.
%                     1 if the intensity process is already known and so its
%                     mean and standard deviation.
%      paths     -    s paths of dimension p*m+1 to calculate a discrete
%                     estimation of the intensity process and then a point
%                     estimator of the intensity mean. s>=30 if the intensity
%                     cannot be assumed as gaussian.
%      timepoints   - p+1 vector (t_0,t_1,...,t_p) of equally or unequally
%                     spaced time points. If they are unequally spaced, at
%                     least delays between observations must be the same
%                     for all m subpaths in which all the s paths will be
%                     splitted.
%      mean_sd      - 2*(p+1) vector. 
%                     The first row are the values of the estimated mean of
%                     the intensity process at each time point (type=0) or the 
%                     values of the actual mean of the intensity process at
%                     each time point (type=1).
%                     The second row are the values of the estimated
%                     standard deviation of the intensity process at each
%                     time point (type=0) or the values of the actual
%                     standard deviation of the intensity process at each
%                     time point (type=1).
%      alpha        - significance level (default = 0.05). 
%  Output:
%      The function returns 1 if data are coherent with the proposal,
%      estimated or real, process and 0 otherwise.
%
%  One example of application of this test with real counting data can be 
%  found in [2], it is applied to data of emitted particles of two different 
%  isotopes. For each new path, a point estimator and an estimated intensity
%  mean and variance (type=0) were taken into account.
%
%
%  Created by (1)P.R. Bouzas (paula@ugr.es) and (2)N. Ruiz-Fuentes
%         (nfuentes@ujaen.es)
%
%         Department of Statistics and Operations Research
%         (1) Universidad de Granada. Campus Universitario de Cartuja. 18071, Granada. Spain
%         (2) Universidad de Jaen. Campus Universitario de Las Lagunillas. 23071, Jaen. Spain  
%
%  Version 1.0. October 2011.
%
%  To cite this file, this would be in an appropriate format:
%  Bouzas, P.R. and Ruiz-Fuentes, N. (2011). Cox process Goodness-of-fit test. 
%  A Matlab file. Free access DIGIBUG
%
%  References:
%  [1] Bouzas, P.R.; Aguilera, A.M.; Valderrama, M.J. and Ruiz-Fuentes, N.,
%              On the structure of the stochastic process of mortgages in
%              Spain, Computation. Stat., 21, 73-89, 2006.
%  [2] Bouzas, P.R.; Ruiz-Fuentes, N.; Matilla, A.;Aguilera, A.M. and
%              Valderrama, M.J., A Cox model for radioactive counting
%              measure: inference on the intensity process, Chemometr.
%              Intell. Lab., 103, 116-121, 2010.
%  [3] Bouzas, P.R.; Aguilera, A.M. and Ruiz-Fuentes, N., Functional
%               estimation of the random rate of a Cox process, Methodol.
%               Comput. Appl. Probab, 2010a, 1-13.
%  [4] Benjamini, Y. and Hochbergk Y., Controling the false rate: a
%               practical and powerful approach to multiple testing, J.
%               Roy. Statist. Soc. Ser. B, 57, 289-300, 1995.

if nargin < 5, 
    alpha = 0.05;  %(default)
end 

if nargin < 4, 
    error('Requires at least four input arguments.'); 
end 

clc
disp('')
disp('----------------------------------------------------------------')
disp('Evaluating if new paths follow the original model')
disp('----------------------------------------------------------------')

T=timepoints;           % observation time points
h=diff(timepoints);     % delays between observation time points
M=paths;                % paths p*m+1 dimensional to calculate a discrete estimation of the intensity process
s=size(M,1);            % number of paths
p=length(h);            
m=(size(M,2)-1)/p;      % number of independent subpaths used for the discrete estimator
Mu=mean_sd(1,:);        % mean of the intensity at each observation point
SD=mean_sd(2,:);        % standard deviation of the intensity at each observation point



%Discrete estimation of the intensity process at the observation points
L=[];
for i=1:s
    L_aux=[];
    for j=1:m
        L_aux(j,:)=diff(M(i,(j-1)*p+1:j*p+1))./h;
    end
    
    L(i,:)=mean(L_aux);
end
Mu_exp=mean(L,1);


%Application of the multiple testing criterion
pv=[];
if type==1
    for j=1:p
        pv(j)=2*(1-normcdf(abs((Mu_exp(j)-Mu(j))/(SD(j)/sqrt(s)))));
    end
elseif type==0
    for j=1:p
        pv(j)=2*(1-tcdf(abs((Mu_exp(j)-Mu(j))/(SD(j)/sqrt(s))),s));
    end
end    
PV=sort(pv');
vale=[]; 
for i=1:p
    if PV(i)>alpha*i/p
       vale(i)=1;
    else
       vale(i)=0;
    end
end  
if sum(vale)==p
   status=1;
   clc
   disp('---------------------------------------------------')
   disp('-----Data are coherent with the proposed model-----')
   disp('----------------------------------------------------')
else
   status=0;
   clc
   disp('-----------------------------------------------------')
   disp('-----Data are incoherent with the proposed model-----')
   disp('-----------------------------------------------------')
end