stem <- function(x,y,pch=16,linecol=1,clinecol=1,...){
# fonction qui permet de dessiner un graphe en bâtons
if (missing(y)){
    y = x
    x = 1:length(x) }
    plot(x,y,pch=pch,...)
       for (i in 1:length(x)){
         lines(c(x[i],x[i]), c(0,y[i]),col=linecol)
        }
   lines(c(x[1]-2,x[length(x)]+2), c(0,0),col=clinecol)
}

plotBinomiale <- function(n, p){
# dessine un histogramme de la loi B(n,p)
  f <- list()
  f$x <- c(0,1)
  nbPoints = length(f$x)
  f$y <- c(1-p, p)
  nf <- loiSomme(f, n)  
  stem(nf$x, nf$y)
  mu = p; sigma2 = p*(1-p);
  t <- seq(min(nf$x), max(nf$x), length.out = 1000);
  lines(t, 1/sqrt(2*pi*n*sigma2)*exp(-(t-n*mu)^2/(2*n*sigma2)), col='red')
}


plotSommeUnif <- function(n){
# dessine la densite de X1+...+Xn, où les Xi sont iid U[0,1]
  f <- list()
  nbPoints <- 101
  f$x <- seq(-1, 1, length.out = nbPoints)
  f$y <- (f$x<0.5)*(f$x>-0.5)+0.0 # loi uniforme sur [-1/2, 1/2]
  nf <- loiSomme(f, n)
  plot(nf, type='l')
  
  mu = 0; sigma2 = 1/12;
#  mu = sum(f$x*f$y)*(f$x[2]-f$x[1]); sigma2 = sum(f$x^2*f$y)*(f$x[2]-f$x[1]) - mu^2
  t <- seq(min(nf$x), max(nf$x), length.out = 1000);
  lines(t, 1/sqrt(2*pi*n*sigma2)*exp(-(t-n*mu)^2/(2*n*sigma2)), col='red')
}

plotGamma <- function(n, l){
# dessine la densite de X1+...+Xn, où les Xi sont iid Exp(l)
  nbPoints <- 201
  nf <- list() 
  nf$x <- seq(n/l*0.1, n/l * 4, length.out = nbPoints)
  nf$y <- dgamma(nf$x, shape=n, rate = l, log = FALSE)
  plot(nf, type='l')
  
  
  
  mu = 1/l; sigma2 = 1/l^2;
  lines(nf$x, 1/sqrt(2*pi*n*sigma2)*exp(-(nf$x-n*mu)^2/(2*n*sigma2)), col='red')  
}



loiSomme <- function(f, n){
# prend en entree la densite (ou proba discrete) f$x,f$y et sort la convolee n fois  
  nf <- f;
  if (n>1) {for (j in seq(2, n)){
    nf$x <- seq(min(nf$x)+min(f$x), max(nf$x)+max(f$x), length.out = length(nf$x) + length(f$x)-1)
    nf$y = convolve(nf$y, rev(f$y), type='open')
    nf$y = nf$y / ((nf$x[2]-nf$x[1])*sum(nf$y))
  }}
  nf
}


#plotBinomiale(50, 0.25)
#plotSommeUnif(4)
#plotGamma(10, 1/2)

tclAPI = function (){
  require(tcltk)
 
#  bernoulli <- tclVar(1)
#  uniform <- tclVar(0)
#  exponential <- tclVar(0)
  
  distrib <- tclVar("bernoulli")
  
  n <- tclVar(5) # nombre de v.a. à ajouter
  pmax <- 100; p <- tclVar(pmax/2) # paramètre de la Bernoulli
  lresol <- 10; l <- tclVar(lresol) # paramètre de l'exponentielle

  base <- tktoplevel() 
  tkwm.title(base,'Central Limit Theorem')
  tt <- tkframe(base, borderwidth=2)
  sliderFrame <- tkframe(tt, borderwidth=2, relief='groove')

  nFrame <- tkframe(sliderFrame, borderwidth=2)
  nLabel <- tklabel(nFrame, text="n ", width=10)

  pFrame <- tkframe(sliderFrame, borderwidth=2)
  pLabel <- tklabel(pFrame, text="p ", width=6)
  pValueLabel <- tklabel(pFrame, text="", width=4)
  
  lFrame <- tkframe(sliderFrame, borderwidth=2)
  lLabel <- tklabel(lFrame, text="l ", width=6)
  lValueLabel <- tklabel(lFrame, text="", width=4)

  actionPerformed = function(...){
    if (tclvalue(distrib)=='bernoulli') {plotBinomiale(as.integer(tclvalue(n)),  as.integer(tclvalue(p))/pmax);}
    else if (tclvalue(distrib)=='uniform') {plotSommeUnif(as.integer(tclvalue(n)));}
    else if (tclvalue(distrib)=='exponential') plotGamma(as.integer(tclvalue(n)), as.integer(tclvalue(l))/lresol);
    tkconfigure(pValueLabel, text=as.character(as.integer(tclvalue(p))/pmax));
    tkconfigure(lValueLabel, text=as.character(as.integer(tclvalue(l))/lresol));
  }

  nSlider <- tkscale(nFrame, from=1, to=25, length=600, showvalue=T, variable=n, resolution=1, orient='horizontal', command=actionPerformed)
  tkpack(nFrame, nLabel, nSlider, side='left');
  pSlider <- tkscale(pFrame, from=0, to=pmax-1, length=600, showvalue=F, variable=p, resolution=1, orient='horizontal', command=actionPerformed)
  tkpack(pFrame, pLabel, pValueLabel, pSlider, side='left');
  lSlider <- tkscale(lFrame, from=1, to=10*lresol, length=600, showvalue=F, variable=l, resolution=1, orient='horizontal', command=actionPerformed)
  tkpack(lFrame, lLabel, lValueLabel, lSlider, side='left');
  tkpack(sliderFrame, nFrame, pFrame, lFrame, side='top')

  rbFrame <- tkframe(tt, borderwidth=2,relief='groove')
  rbBernoulli <- tkradiobutton(rbFrame, command=actionPerformed, text='Bernoulli', variable=distrib)
  rbUniform <- tkradiobutton(rbFrame, command=actionPerformed, text='Uniform', variable=distrib)
  rbExponential <- tkradiobutton(rbFrame, command=actionPerformed, text='Exponential', variable=distrib)
  tkconfigure(rbBernoulli, variable=distrib, value='bernoulli'); 
  tkconfigure(rbUniform, variable=distrib, value='uniform');
  tkconfigure(rbExponential, variable=distrib, value='exponential');
  tkpack(rbFrame, rbBernoulli, rbUniform, rbExponential, side='left');
  
  tkpack(tt, rbFrame, sliderFrame, side='top');
  actionPerformed()
}


tclAPI();