%
% CpTruncNorm Matlab file, v.1.0.2011
%
%  CpTruncNorm is inspired in the outcomes in Bouzas et al. [1], it
%  calculates the most probable value of a Cox process with mean a
%  truncated Gaussian distribution. It was proved (Bouzas et al.[1]) that
%  this value, the mode of the DSPP at each time point, is within a bounded
%  interval and can be calculated.
%
% Syntax: function [Mo,pr,cdf]=TruncNorm(A,B,M,S,l)
%
%  Inputs:
%      A         - Lower truncation point of the truncated Gaussian distribution.
%      B         - Upper truncation point of the truncated Gaussian distribution.
%      M         - Mean of the truncated Gaussian distribution.
%      S         - Standard deviation of the truncated Gaussian distribution.
%      l         - The distribution function of the specified truncated Gaussian 
%                  distribution will be evaluated at this point.Point at which 
%
%  Output:
%      It returns Mo, that is the most probable value of the Cox process in
%      an instant of time with mean a truncated Gaussian distribution with
%      the specified parameters at this instant. It returns pr, it is the
%      probability of ocurrence of the previous most probable value. cdf is
%      the value of the distribution function of the truncated Gaussian
%      random variable at l, this is an additional task solved by the
%      function.
%
%  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. November 2011.
%
%  To cite this file, this would be in an appropriate format:
%  Bouzas, P.R. and Ruiz-Fuentes, N. (2011). CpTruncNorm Matlab file. Free access DIGIBUG
%
%  References:
%  [1] Bouzas, P.R., A.M. Aguilera and Valderrama, M.J., Forecasting a class of 
%              doubly stochastic Poisson process, Statistical Papers, 43, 
%              507--523, 2002.
function [Mo,pr,cdf]=CpTruncNorm(A,B,M,S,l)
clc
disp('                                              ')
syms s L u
k1=(S*int(exp(-s^2/2),(A-M)/S,(B-M)/S))^-1;
f=inline((S*int(exp(-s^2/2),(A-M)/S,(B-M)/S))^-1*exp((-1/2)*((L-M)/S)^2));
x=[A:(B-A)/100:B];
plot(x,f(x));title('Density function');

f1=inline(int(exp(-s^2/2),(A-M)/S,(u-M)/S));
cdf=f1(l)/f1(B);
disp('                                              ')
disp('----------------------------------------------')
disp('The distribution function at the given point is '),disp([cdf]),

k2=(1/sqrt(2*pi))*f1(B);
P=[];
syms s
for n=floor(A):floor(B)
    p=inline(int(s^n*exp(-s)*exp(-0.5*(s-M)^2/S^2),A,B));
    P=[P;(1/factorial(n))*k2^-1*(1/sqrt(S^2*2*pi))*p(n)];
end
[pr,i]=max(P);
Mo=floor(A)+i;
disp('The maximun probability integer is '),disp(floor(A)+i),
disp('With probability                   '),disp(pr);