#  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 Bernstein polynomials    ==== 

#computing a vector of lenght n+1 with values of Bernstein polynomials at x
bernstein_poly <- function( n, x ) .Call("bernstein_poly_vec", n, x) 

#computing quantiles using B-polynomials for the empirical quantile function
bpq_quantile=function(sample, probs, k=100)
{
 q=probs;
 nonsm_q=quantile(sample,c(1:k)/k,type=1);
 for(i in (1:length(probs)))
 {
    q[i]=sum(bernstein_poly ( k-1, probs[i] )*nonsm_q);
 }
 return(q)
}

#computing values of the emprical distribution function at specified points
edf_at=function(sample,grid)
{
 vals=grid;
 m=length(sample);
 for( i in (1:length(grid)) )
 {
    vals[i]=sum(sample<grid[i])/m;
 }
 return(vals)
}

#computing quantiles using B-polynomials for the empirical distribution function
bpf_quantile=function(sample, probs, k=100)
{
 grid01 = c(1:k)/(k+1);
 grid = tan((grid01-0.5)*pi);
 a=mean(sample);
 s=sd(sample);
 edf = edf_at((sample-a)/s, grid);
 cdf = bpq_quantile(edf, grid01, k);
 q=probs;
 for(i in (1:length(probs)))
 {
    kk=max(which(cdf<probs[i]));
    gam=(probs[i]-cdf[kk])/(cdf[kk+1]-cdf[kk]);
    t=grid01[kk]+gam*(grid01[kk+1]-grid01[kk]);
    q[i]=a+s*tan((t-0.5)*pi);
 }
 return(q)
}

#R-computing a vector of lenght n+1 with values of Bernstein polynomials at x
bernstein_poly_slow=function( n, x )
{
  n=abs(n);
  if(n == 0)
  {
    bern = c(1.0);
  }
  else
  {  
    bern=rep(0,n+1);
    bern[1] = 1.0 - x;
    bern[2] = x;
    if(n>1)
    {
    for(i in c(2:n))
    {
      bern[i+1] = x * bern[i];
      for( j in c((i-1):1) )
      {
        bern[j+1] = x * bern[j] + ( 1.0 - x ) * bern[j+1];
      }
      bern[1] = ( 1.0 - x ) * bern[1];
    }
    }
  }
  return(bern);
}