require(MASS);

epanechnikov <- function (t) {3/4*(1-t^2)*(t>-1)*(t<1)}

rePlot = function(typeGraph, ...){#y, typeSource, lambda, plotRef) {
  if (typeGraph=='Histo') {plotHist(...);}
  else if (typeGraph=='Kernel') {plotKernel(...);}
  else {plotFR(...);}
}

plotHist = function(y, lambda, typeSource, plotRef, ...){
  truehist(y, xlim=c(-4, 4), ylim=c(0,0.6), nbins=max(1, floor(6/lambda)), xlab='empirical measure', ylab='histogram', col='blue');#nbins='FD');
  if (length(y)<501) points(y, array(-0.01, length(y)), pch=19, col='cyan', lwd=5)
  if (plotRef){
    xx<-seq(-5, 5, 0.01);
    if (typeSource=='Gaussian'){lines(xx, dnorm(xx), col='red', lwd=2);}
    else if (typeSource=='Mixt2') {lines(xx, dnorm(2*(xx-2))+dnorm(2*(xx+2)), col='red', lwd=2);}
    else if (typeSource=='Mixt3') {lines(xx, 2/3*dnorm(2*xx)+2/3*dnorm(2*(xx-2))+2/3*dnorm(2*(xx+2)), col='red', lwd=2);}
  }
}

plotKernel = function(y, lambda, typeSource, plotRef, kernelType){
  k<-length(y);
  tt <- seq(-5, 5, 0.01);
  fh <- tt-tt; 
  plot(0,0,type='n', xlim=c(-4, 4), ylim=c(0,0.6), xlab='empirical measure', ylab='density');
  if (k<501) lines(y[1:k], array(0, k), pch=19, col='cyan', lwd=5);
  for (j in 1:k){
      if (kernelType=='Epanechnikov'){
      contrib <- 1/(k*lambda) * epanechnikov((tt-y[j])/lambda);
    }
      else{ # Gaussian Kernel
      contrib <- 1/(k*lambda) * dnorm((tt-y[j])/lambda);
    }
    fh <- fh + contrib;
    if (k<101) lines(tt, contrib, col='green'); #  lines(x[j]*c(1,1), c(0, max(contrib)), type='l', col='black', lwd=2);
  }
  lines(tt, fh, type='l', lwd=3, col='blue');
  if (plotRef){
    if (typeSource=='Gaussian'){lines(tt, dnorm(tt), col='red', lwd=2);}
    else if (typeSource=='Mixt2') {lines(tt, dnorm(2*(tt-2))+dnorm(2*(tt+2)), col='red', lwd=2);}
    else if (typeSource=='Mixt3') {lines(tt, 2/3*dnorm(2*(tt-2))+2/3*dnorm(2*tt)+2/3*dnorm(2*(tt+2)), col='red', lwd=2);}
  }
}

plotFR = function(y, lambda, typeSource, plotRef,...){
  plot(c(-4, sort(y), 4), c((0:length(y))/length(y), 1), type='s', col='blue', xlim=c(-4, 4), ylim=c(0,1.01), lwd=2, xlab='empirical measure', ylab='cumulative distribution function');#nbins='FD');
  if (length(y)<501) points(y, array(-0.01,length(y)), pch=19, col='cyan', lwd=5)
  if (plotRef){
    xx<-seq(-5, 5, 0.01);
    if (typeSource=='Gaussian'){lines(xx, pnorm(xx), col='red', lwd=2);}
    else if (typeSource=='Mixt2') {lines(xx, 1/2*pnorm(2*(xx-2))+1/2*pnorm(2*(xx+2)), col='red', lwd=2);}
    else if (typeSource=='Mixt3') {lines(xx, 1/3*pnorm(2*(xx-2))+1/3*pnorm(2*xx)+1/3*pnorm(2*(xx+2)), col='red', lwd=2);}
  }
}