hyp.test<-function(n=5,sd=1,nreps=100) {
  ci<-matrix(numeric(40),ncol=2,nrow=nreps);
  pval<-numeric(nreps);
  for(i in 1:nreps) {
    a<-rnorm(n,10,sd);
    b<-rnorm(n,10,sd);
    tt<-t.test(a,b,var.equal=TRUE);
    ci[i,]<-tt$conf.int[1:2];
    pval[i]<-tt$p.value;
  }

  par(mfrow=c(2,1));
  ci.lim<-range(ci);
  x<-1:nreps;
  par(mar=c(0,4,1,1)+0.1);
  plot(x,rep(0,nreps),type="l",lty=2,col="red",xlab="Tests de t",ylab="Différence de moyenne",ylim=ci.lim,xaxt="n");
  for(i in 1:nreps) {
    segments(x[i],ci[i,1],x[i],ci[i,2],col="blue")
  }
  par(mar=c(4,4,1,1)+0.1);
  plot(x,pval,ylim=c(0,1),xlab="Tests de t",ylab="P-value",pch="x");
  lines(x,rep(0.05,nreps),col="red");

#return(list(ci=ci,p=pval));
}

simulH1<-function(n=5,ne=0,sd=1,unil=FALSE,plot=TRUE,power=TRUE) {
  if(n<2) {
    n<-2;
    cat("Nombre de données a été fixé à 2\n");
  }

  if(ne<0) {
    ne<-abs(ne);
    cat("Différence entre H0 et H1 ne peut pas être négative, valeur absolue prise\n");
  }
  df<-(n-1)*2;
  ncp<-(ne/sd)*sqrt(n);
  b<-rnorm(n,0,sd);
  a<-rnorm(n,ne,sd);
  if(plot) {
    vec<-seq(-10,ncp+10,0.001);
    density.H0<-dt(vec[vec<10],df);
    if(ne>0) {
      lim.H1<-ncp-4*sd;
      density.H1<-dt(vec[vec>=lim.H1],df,ncp);
      my.lim<-c(
        min(vec[vec<10][min(which(density.H0>0.001))],vec[vec>=lim.H1][min(which(density.H1>0.001))]),
        max(vec[vec<10][max(which(density.H0>0.001))],vec[vec>=lim.H1][max(which(density.H1>0.001))])
      );
    }
    else {
      my.lim<-vec[range(which(density.H0>0.001))];
    }
    par(mar=c(4,4,1,1)+0.1);
    plot(vec[vec<10],density.H0,type="l",col="red",xlab="Statistique de T",ylab="densité", xlim=my.lim);
    if(ne>0) {
      lines(vec[vec>=lim.H1],density.H1,type="l",col="blue");
    }
    if(unil) {
      perc<-qt(0.95,df);
      if(power) {
        power.t<-1-pt(perc,df,ncp);
      }
      else {
        power.t<-1.0;
      }
      tt<-t.test(a,b,var.equal=TRUE,alternative="greater");
      wt<-wilcox.test(a,b,alternative="greater");
      vec.1<-seq(perc,my.lim[2],0.001);
      if(ne>0) {
        lines(vec.1,dt(vec.1,df,ncp),type="h",col="blue");
      }
      lines(vec.1,dt(vec.1,df),type="h",col="red");
    }
    else {
      tt<-t.test(a,b,var.equal=TRUE);
      wt<-wilcox.test(a,b);
      perc<-qt(0.975,df);
      if(power) {
        power.t<-1-pt(perc,df,ncp);
      }
      else {
        power.t<-1.0;
      }
      vec.1<-seq(perc,my.lim[2],0.001);
      if(ne>0) {
        lines(vec.1,dt(vec.1,df,ncp),type="h",col="blue");
      }
      lines(vec.1,dt(vec.1,df),type="h",col="red");
      perc<-qt(0.025,df);
      vec.1<-seq(my.lim[1],perc,0.001);
      lines(vec.1,dt(vec.1,df),type="h",col="red");
    }
    pt<-tt$p.value;
    t<-tt$statistic;
    ci.t<-tt$conf.int[1:2]
    pw<-wt$p.value;
    w<-wt$statistic;
    abline(h=0);
    abline(v=t,col="green");
    if(ne>0) {
      height<-max(density.H1)/2;
      x<-(my.lim[2]-my.lim[1])*0.3/13; #the ideal size of the white box behind H1, changing with xlim
      rect(ncp-x,height-0.01,ncp+x,height+0.01,col="white",border="white");
      rect(-x,height-0.01,x,height+0.01,col="white",border="white");
      text(ncp,height,labels="H1",col="blue");
      text(0,height,labels="H0",col="red");
    }
    else {
      height<-max(density.H0)/2;
      x<-(my.lim[2]-my.lim[1])*0.3/13; #the ideal size of the white box behind H1, changing with xlim
      rect(-x,height-0.01,x,height+0.01,col="white",border="white");
      text(0,height,labels="H0",col="red");
    }
    legend(my.lim[1],max(density.H0),legend="T obs.",pch="_",col="green");
  }
  else {
    if(unil) {
      tt<-t.test(a,b,var.equal=TRUE,alternative="greater");
      wt<-wilcox.test(a,b,alternative="greater");
      if(power) {
        perc<-qt(0.95,df);
        power.t<-1-pt(perc,df,ncp);
      }
      else {
        power.t<-1.0;
      }
    }
    else {
      tt<-t.test(a,b,var.equal=TRUE);
      wt<-wilcox.test(a,b);
      if(power) {
        perc<-qt(0.975,df);
        power.t<-1-pt(perc,df,ncp);
      }
      else {
        power.t<-1.0;
      }
    }
    pt<-tt$p.value;
    t<-tt$statistic;
    ci.t<-tt$conf.int[1:2]
    pw<-wt$p.value;
    w<-wt$statistic;
  }
  if(power) {
    power.w<-power.wilcox.test(ne=ne,n=n,sd=sd);
  }
  else {
    power.w<-1.0;
  }
#res<-matrix(c(t,pt,power.t,w,pw,power.w),ncol=3,nrow=2,byrow=TRUE);
  res<-matrix(c(t,pt,ci.t),ncol=4,nrow=1,byrow=TRUE);
  rownames(res)<-c("Test de T");
  colnames(res)<-c("Statistiques","P-value","Int. conf. inf.", "Int. conf. sup.");
  return(list(res=res));
}

pvalue<-function(nreps=100,p=0.05,n=5,sd=1,unil=FALSE) {
  mat<-matrix(numeric(2*nreps),nrow=nreps,ncol=2);
  for(i in 1:nreps) {
    mat[i,]<-simulH1(n,0,sd,unil,FALSE,FALSE)$res[,2];
  }
  res<-length(which(mat[,1]<p));
#cat("Sur ",nreps," tests (n=",n,", sd=",sd,"), ",res," sont significatifs.\n",sep="");
  return(res);
}

power.wilcox.test<-function(ne=1,n=5,sd=1,unil=FALSE) {
  wilcox.pval<-function(ne=1,n=5,sd=1,unil=FALSE) {
    a<-rnorm(n,0,sd);
    b<-rnorm(n,ne,sd);
    if(unil) {
      wilcox.test(a,b,alternative="greater")$p.value;
    }
    else {
      wilcox.test(a,b)$p.value;
    }
  }

  pvals<-replicate(1000,wilcox.pval(ne,n,sd));
  return(sum(pvals<0.05)/1000);
}

independent<-function(n=5,ne=0.5,sd=1,unil=FALSE) {
  a<-rnorm(n/2,0,sd);
  a<-c(a,a);
  b<-rnorm(n/2,ne,sd);
  b<-c(b,b);
  if(unil) {
    tt<-t.test(a,b,var.equal=TRUE,alternative="greater");
    wt<-wilcox.test(a,b,alternative="greater");
  }
  else {
    tt<-t.test(a,b,var.equal=TRUE);
    wt<-wilcox.test(a,b);
  }
  pt<-tt$p.value;
  t<-tt$statistic;
  ci.t<-tt$conf.int[1:2]
  pw<-wt$p.value;
  w<-wt$statistic;
  res<-c(t,pt,w,pw);
  names(res)<-c("T.stat", "T.pval", "Wilcox.stat", "Wilcox.pval");
  return(res);
}

normality<-function() {
}