#  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 2 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.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/

### copyright (C) 2011 A. Pepelyshev
### This version distributed under GPL (version 2 or later)


#==== nonparametric estimation using orthogonal series expansion  ==== 

#estimating the density and cumulative distribution function
ose_density=function(sam,type=1,N=10)
{
# type 1,2,3,4 = Kronmal-Tarter 1,2,3,4
# type -1 = modified Kronmal-Tarter (include nonzero coefficients)
# type -2 = Diggle-Hall
KT=type;
m=length(sam);
if(type<0)
{
    KT=-type;
}
if(KT>0) # Kronmal-Tarter (t)
{
    N=min(100,round(m/3));
}
if(type==-2) # Diggle-Hall
{
       N=310;
}

a=min(sam)-0.5*sd(sam);
b=max(sam)+0.5*sd(sam);
xab=(sam-a)/(b-a);

cc=rep(1,N); d2=cc;
q=sqrt(2);
for(i in (2:N))
{
    cc[i]=mean(cos(pi*(i-1)*xab))*q;
    d2[i]=mean( (cos(pi*(i-1)*xab)*q)^2 );
}
Jset=c(2:N);
if(type>0) # Kronmal-Tarter (t)
{
    isf=0;
    for(i in (4:(N-5)) )
    {
        if((isf==0) && (sum(cc[(i+1):(i+KT)]^2<d2[(i+1):(i+KT)]*2/(m+1))==KT))
        {
            Jset=c(2:i);
            isf=1;
        }
    }
}
if(type==-1) # modified Kronmal-Tarter
{
    Jset=c(0);
    for(i in (2:N))
    {
        if(cc[i]^2>d2[i]*2/(m+1))
        {
            Jset=c(Jset,i);
        }
    }
    Jset=Jset[-1];
}
if(type==-2) # Diggle-Hall
{
    k=floor((N/3)^(2/3));
    psi=rep(1e+30,k);
    for(i in (4:k))
    {
        ii=floor(3*i^(3/2));
        psi[i]=i/(m*pi)+sum( max(0*c(i:ii),cc[i:ii]^2-d2[i:ii]/m) );
    }
    ii=which.min(psi);
    Jset=c(2:ii);
}


K=1000;
tm=c(0:K)/K;
appr=rep(0,length(tm))+cc[1];
for(i in Jset)
{
    appr=appr+cc[i]*q*cos(pi*(i-1)*tm);
}

isOK=1;
while(isOK)
{
    ap=appr-0.001;
    ap[ap<0]=0;
    integ=simp1(tm,ap);
    if(integ$i>1)
    {
        appr=ap;
    }
    else
    {
        isOK=0;
    }
}
integ=simp1(tm,appr);
s=integ$i;
ecdf=integ$cdf;

return(list(tm=(a+tm*(b-a)),cdf=ecdf,pdf=appr))
}

#computing quantiles
ose_quantile=function(sam,probs,type=1,N=10)
{
ose=ose_density(sam,type,N);
q=probs;
for(i in (1:length(probs)) )
{
    k=max(which(ose$cdf<probs[i]));
    alp=(probs[i]-ose$cdf[k])/(ose$cdf[k+1]-ose$cdf[k]);
    q[i]=ose$tm[k]+alp*(ose$tm[k+1]-ose$tm[k]);
}
return(q);
}


simp1=function(x,y)
{
#SIMP  Simpson's method numerical integration.
# Sum of parabolic arcs computed with vector-matrix multiply.
m = length(y);
c1 = rep(1,m-1); c1[seq(2,m-1,2)] = 2; 
c2 = rep(1,m-1); c2[seq(1,m-1,2)] = 2; 
csum = c(0,cumsum(diff(x)*(c1*y[1:m-1] + c2*y[2:m])/3));
z = csum[m];
return(list(i=z,cdf=csum))
}