#  (gH)   --  lois.pl  ;  TimeStamp (unix) : 09 Janvier 01 11:45

sub syntaxe {

  $ligSyntax  = " [ b p | B n p | G p t | H n k s | P m t |  Q s p t \n" ;
  $ligSyntax .= (" "x18)."      | T a   | U a b | S m t   | Z ] T " ;

  print "\n" ;
  print " syntaxe  : lois $ligSyntax\n\n" ;
  print " exemples : \n" ;
  print "   lois  b   0.2       # loi de bernoulli     p=0.2\n" ;
  print "   lois  B   5 0.1     # loi binomiale        n=5, p=0.1\n" ;
  print "   lois  G   0.8 9     # loi géométrique      p=0.8 limitée à 9 termes\n" ;
  print "   lois  H   10 5 3    # loi hypergéométrique N=10, K=5, n=3\n" ;
  print "   lois  N   4 0.9 11  # loi de pascal        s=4, p=0.9 limitée à 11 termes\n" ;
  print "   lois  P   1.2 9     # loi de poisson       m=1.2 limitée à 9 termes\n" ;
  print "   lois  S   0.3 9     # loi \"sans-nom\"       p=0.3 limitée à 9 termes\n" ;
  print "   lois  T   4         # loi triangulaire     a=4\n" ;
  print "   lois  U   2   7     # loi uniforme         a=2, b=7\n" ;
  print "   lois  Z             # loi de zipf-benford           \n" ;
  print "\n T est la taille de l'échantillon correspondant (100 par défaut)." ;
  print "\n" ;

} ; # fin sub syntaxe

##########################################################################
##                                                                      ##
##   programme principal                                                ##
##                                                                      ##
##########################################################################

print "\n" ;
print " (gH) ; lois.pl v1.2 : lois statistiques discrètes usuelles \n" ;


if ($#ARGV==-1) { &syntaxe() ; exit 0 ; }  ;

$nomLoi{"b"} = "de bernoulli" ;
$nomLoi{"B"} = "binomiale (compte-avec)" ;
$nomLoi{"G"} = "géométrique (bernoulli négative)" ;
$nomLoi{"H"} = "hypergéométrique (compte-sans)" ;
$nomLoi{"N"} = "de pascal (binomiale négative ou première-avec)" ;
$nomLoi{"P"} = "de poisson" ;
$nomLoi{"S"} = "sans-nom (première-sans)" ;
$nomLoi{"U"} = "uniforme discrète" ;
$nomLoi{"T"} = "triangulaire" ;
$nomLoi{"Z"} = "de zipf-benford" ;

if ($ARGV[0] eq "b") { &bernoulli($ARGV[1],$ARGV[2])                           ; exit 0 ; }  ;
if ($ARGV[0] eq "B") { &binomiale($ARGV[1],$ARGV[2],$ARGV[3])                  ; exit 0 ; }  ;
if ($ARGV[0] eq "G") { &geometrique($ARGV[1],$ARGV[2],$ARGV[3])                ; exit 0 ; }  ;
if ($ARGV[0] eq "H") { &hypergeometrique($ARGV[1],$ARGV[2],$ARGV[3],$ARGV[4])  ; exit 0 ; }  ;
if ($ARGV[0] eq "N") { &pascal($ARGV[1],$ARGV[2],$ARGV[3],$ARGV[4])            ; exit 0 ; }  ;
if ($ARGV[0] eq "P") { &poisson($ARGV[1],$ARGV[2],$ARGV[3])                    ; exit 0 ; }  ;
#if ($ARGV[0] eq "S") { &sans($ARGV[1],$ARGV[2])                      ; exit 0 ; }  ;
if ($ARGV[0] eq "T") { &triangulaire($ARGV[1],$ARGV[2])                        ; exit 0 ; }  ;
if ($ARGV[0] eq "U") { &uniforme($ARGV[1],$ARGV[2],$ARGV[3])                   ; exit 0 ; }  ;
if ($ARGV[0] eq "Z") { &zipfbenford($ARGV[1])                                  ; exit 0 ; }  ;

print "\n loi $ARGV[0] : pas encore implémentée\n         ou faute de frappe sur l'initiale de la loi !\n" ;
exit 0 ;

#########################################################################

sub log10 { return log($_[0])/log(10) ; } ;

#########################################################################

sub cnk { # coefficients du binome

  my ($n,$p) = ($_[0],$_[1]) ;
  if ($p>$n) { return 0 ; } ;
  if ($n<0)  { return 0 ; } ;
  if ($p<0)  { return 0 ; } ;

  $k = $p ;
  if ($k>$n-$k) { $k = $n-$k ; } ;
  if    ($k==0) { $res = 1  ; }
  elsif ($k==1) { $res = $n ; }
  elsif ($k==2) { $res = $n*($n-1)/2 ; }
  else  {
    $i    = 0.0 ;
    $frac = 1.0 ;
    while ($i<$p) {
      $i++ ;
      $frac = $frac * ($n+1-$i)/$i ;
    } ; # fin tant que
    $res = $frac ;
  } ; # fin de si
  ## print "c($n,$p) = $res \n" ;
  return $res ;
} ; # fin sub cnk

#########################################################################

sub bernoulli {

  $proba     = $_[0] ;
  $taille    = $_[1] ;
  $nbVal     = 2 ;
  $paramLoi  = "p=$proba" ;
  ($tabVal[1],$tabVal[2]) = (0,1) ;
  ($tabPro[1],$tabPro[2]) = (1-$proba,$proba) ;

  &decritLoi("b") ;

} ; # fin sub bernoulli

#########################################################################

sub binomiale {

  my $n      = $_[0]  ;
  $nbVal     = 1 + $n ;
  $proba     = $_[1]*1.0  ;
  $taille    = $_[2] ;
  $qroba     = 1 - $proba ;
  $paramLoi  = "n=$n, p=$proba" ;
  my $i = 0 ;
  while ($i<=$n) {
    $tabVal[1+$i] = $i ;
    $tabPro[1+$i] = &cnk($n,$i)*($proba**$i)*($qroba**($n-$i)) ;
    $i++;
  } ; # fin tant que

  &decritLoi("B") ;

} ; # fin sub binomiale

#########################################################################

sub uniforme {

  $a         = $_[0]  ;
  $b         = $_[1]  ;
  $taille    = $_[2] ;
  $nbVal     = 1 + $b - $a ;
  $proba     = 1.0/$nbVal  ;
  $paramLoi  = "a=$a, b=$b" ;
  my $i = 1 ;
  while ($i<=$nbVal) {
    $tabVal[$i] = $a + ($i-1) ;
    $tabPro[$i] = $proba ;
    $i++;
  } ; # fin tant que

  &decritLoi("U") ;

} ; # fin sub uniforme

#########################################################################

sub poisson {

  $m         = $_[0]  ;
  $nbVal     = $_[1]  ;
  $taille    = $_[2] ;
  $puism     = 1.0 ;
  $factk     = 1.0 ;
  $paramLoi  = "m=$m (limitée à $nbVal termes)" ;
  my $i   = 1 ;
  my $kst = exp(-$m) ;
  while ($i<=$nbVal) {
    $tabVal[$i] = $i-1 ;
    $tabPro[$i] = $kst*$puism/$factk ;
    $factk = $factk * $i ;
    $i++;
    $puism = $puism * $m ;
  } ; # fin tant que

  &decritLoi("P") ;

} ; # fin sub poisson

#########################################################################

sub geometrique {

  $p         = $_[0]   ;
  $q         = 1 -$p ; ;
  $nbVal     = $_[1]   ;
  $taille    = $_[2] ;
  $paramLoi  = "p=$p (limitée à $nbVal termes)" ;
  my $i      = 1 ;
  while ($i<=$nbVal) {
    $tabVal[$i] = $i ;
    $tabPro[$i] = $p*($q**($i-1)) ;
    $i++;
  } ; # fin tant que

  &decritLoi("G") ;

} ; # fin sub geometrique

#########################################################################

sub pascal {

  $s         = $_[0]   ;
  $p         = $_[1]   ;
  $q         = 1.0 - $p ; ;
  $nbVal     = $_[2]   ;
  $taille    = $_[3] ;
  $paramLoi  = "s=$s, p=$p (limitée à $nbVal termes)" ;
  my $i   = 1 ;
  while ($i<=$nbVal) {
    $tabVal[$i] = $i ;
    if ($i<$s) { $tabPro[$i] = 0 ; } else {
       $j = $i-$s ;
       $x = &cnk($s+$j-1,$j) ;
       $y = $p**$s ;
       $z = $q**$j ;
       $lap = $x * $y * $z ;
       $tabPro[$i] = $lap ;
    } ; # fin de si
    $i++;
  } ; # fin tant que

  &decritLoi("N") ;

} ; # fin sub pascal

#########################################################################

sub hypergeometrique {

  $n1        = $_[0]   ;
  $n2        = $_[1]   ;
  $n3        = $_[2]   ;
  $taille    = $_[3] ;
  $paramLoi  = "N=$n1, K=$n2, n=$n3" ;
  $a         = $n3 - ($n1-$n2) ;
  if ($a<0)   { $a = 0 } ;
  $b         = $n3 ;
  if ($b>$n2) { $b = $n2 } ;
  $nbVal = $b - $a + 1 ;
  my $i      = 1 ;
  while ($i<=$nbVal) {
    $j          = $a + $i - 1 ;
    $tabVal[$i] = $j ;
    $tabPro[$i] = &cnk($n2,$j)*&cnk($n1-$n2,$n3-$j)/&cnk($n1,$n3) ;
    $i++;
  } ; # fin tant que

  &decritLoi("H") ;

} ; # fin sub hypergeometrique

#########################################################################

sub triangulaire {

  $a         = $_[0]   ;
  $taille    = $_[1] ;
  $paramLoi  = "a=$a" ;
  $nbVal = 2*$a + 1 ;
  my $i      = 1 ;
  while ($i<=$nbVal) {
    $j          = -$a + $i - 1 ;
    $tabVal[$i] = $j ;
    $tabPro[$i] = ($a+1-abs($j))/(($a+1)*($a+1));
    $i++;
  } ; # fin tant que

  &decritLoi("T") ;

} ; # fin sub triangulaire

#########################################################################

sub zipfbenford {

  $paramLoi  = " (aucun)" ;
  $taille    = $_[0] ;
  $i = 1 ;
  $nbVal = 9 ;
  while ( $i <= $nbVal) {
    $tabVal[$i] = $i ;
    $tabPro[$i] = log10(1+$i) - log10($i) ;
    $i++;
  } ; # fin tant que

  &decritLoi("Z") ;

} ; # fin sub zipfbenford

#########################################################################

sub decritLoi {

  print "\nloi ".$nomLoi{$_[0]}.", paramètre(s) : $paramLoi\n" ;
  if ($taille eq "") { $taille = 100 } ;
  print "   valeurs        probabilité   cumul      effectif (taille =$taille)\n" ;
  my $iv  = 1 ; # indice de valeur
  my $cpr = 0.0 ; # cumul des probabilités
  my $spv = 0.0 ; # somme pondérée des valeurs
  my $spc = 0.0 ; # somme pondérée des carrés des valeurs
  my $sde = 0.0 ; # somme des effectifs

  $lsor = "lois.sor" ;
  open(LSOR,">$lsor") or die("impossible d'écrire dans $lsor\n") ;

  while ($iv<=$nbVal) {
    $vc   = $tabVal[$iv] ;
    $vp   = $tabPro[$iv] ;
    $vf   = $vp*$taille ;
    $cpr += $vp ;
    $spv += $vc*$vp ;
    $spc += $vc*$vc*$vp ;
    $fmt_v = sprintf("%4d",$vc) ;
    $fmt_p = sprintf("%8.5f",$vp) ;
    $fmt_c = sprintf("%8.5f",$cpr) ;
    $fmt_f = sprintf("%5.0f",$vf) ;
    $sde  += $fmt_f ;
    print "  $fmt_v           $fmt_p      $fmt_c $fmt_f\n" ;
    print LSOR "  $fmt_v           $fmt_p      $fmt_c $fmt_f\n" ;
    $iv++;
  } ; # fin tant que sur $iv
  close(LSOR) ;

  my $moy = $spv ;
  my $var = $spc - $moy*$moy ;
  my $ect = sqrt($var) ;
  if ($moy==0) { $cdv = -1 ; } else { $cdv = sprintf("%4.0f",100.0*$ect/$moy) ; } ;
  $fmt_moy = sprintf("%9.4f",$moy) ;
  $fmt_ect = sprintf("%9.4f",$ect) ;
  print " moyenne calculée  $fmt_moy   écart-type calculé   $fmt_ect  c.d.v. $cdv % \n" ;

  if ($ARGV[0] eq "b") {
     $mot = $proba ;
     $ett = sqrt($proba*(1-$proba)) ;
     $fmt_mot = sprintf("%9.4f",$mot) ;
     $fmt_ett = sprintf("%9.4f",$ett) ;
     print " moyenne théorique $fmt_mot   écart-type théorique $fmt_ett \n" ;
     print " = p                           = racine(p(1-p))\n" ;
  } ; # fin si sur loi bernoulli

  if ($ARGV[0] eq "B") {
     $mot = ($nbVal-1)*$proba ;
     $ett = sqrt(($nbVal-1)*$proba*$qroba) ;
     $fmt_mot = sprintf("%9.4f",$mot) ;
     $fmt_ett = sprintf("%9.4f",$ett) ;
     print " moyenne théorique $fmt_mot   écart-type théorique $fmt_ett \n" ;
     print " = n.p                         = racine(n.p(1-p))\n" ;
  } ; # fin si sur loi binomiale

  if ($ARGV[0] eq "P") {
     $mot = $m ;
     $ett = sqrt($m) ;
     $fmt_mot = sprintf("%9.4f",$mot) ;
     $fmt_ett = sprintf("%9.4f",$ett) ;
     print " moyenne théorique $fmt_mot   écart-type théorique $fmt_ett \n" ;
     print " = m                           = racine(m)\n" ;
  } ; # fin si sur loi poisson

  if ($ARGV[0] eq "G") {
     $mot = 1.0/$p ;
     $ett = sqrt( (1.0-$p)/($p*$p)  ) ;
     $fmt_mot = sprintf("%9.4f",$mot) ;
     $fmt_ett = sprintf("%9.4f",$ett) ;
     print " moyenne théorique $fmt_mot   écart-type théorique $fmt_ett \n" ;
     print " = 1/p                         = racine( (1-p)/p^2 )\n" ;
  } ; # fin si sur loi géométrique

  if ($ARGV[0] eq "N") {
     $mot = $s*1.0/$p ;
     $ett = sqrt( $s*(1.0-$p)/($p*$p)  ) ;
     $fmt_mot = sprintf("%9.4f",$mot) ;
     $fmt_ett = sprintf("%9.4f",$ett) ;
     print " moyenne théorique $fmt_mot   écart-type théorique $fmt_ett \n" ;
     print " = s/p                         = racine( s(1-p)/p^2 )\n" ;
  } ; # fin si sur loi pascal

  if ($ARGV[0] eq "H") {
     $mot = $n3*$n2*1.0/$n1 ;
     $ett = sqrt( ($n1-$n3)*($n1-$n2)*$n3*$n2/( $n1*$n1*($n1-1)) ) ;
     $fmt_mot = sprintf("%9.4f",$mot) ;
     $fmt_ett = sprintf("%9.4f",$ett) ;
     print " moyenne théorique $fmt_mot   écart-type théorique $fmt_ett \n" ;
     print " = nK/N                        = racine( nK(N-n)(N-K) / (n-1)n^2 )\n" ;
  } ; # fin si sur loi pascal

  if ($ARGV[0] eq "U") {
     $mot = ($a+$b)/2 ;
     $ett = sqrt(($b-$a)*($b-$a+2)/12.0) ;
     $fmt_mot = sprintf("%9.4f",$mot) ;
     $fmt_ett = sprintf("%9.4f",$ett) ;
     print " moyenne théorique $fmt_mot   écart-type théorique $fmt_ett \n" ;
     print " = (a+b)/2                     = racine( (b-a)(b-a+2)/12 )\n" ;
  } ; # fin si sur loi uniforme

  if ($ARGV[0] eq "T") {
     $mot = 0 ;
     $ett = sqrt($a*($a+2)/6.0) ;
     $fmt_mot = sprintf("%9.4f",$mot) ;
     $fmt_ett = sprintf("%9.4f",$ett) ;
     print " moyenne théorique $fmt_mot   écart-type théorique $fmt_ett \n" ;
     print "                               = racine( n(n+2)/6)\n" ;
  } ; # fin si sur loi triangulaire

  $f_sde = sprintf("%5.0f",$sde) ;
  print " somme des effectifs arrondis : $f_sde (taille $taille)\n" ;
  print "\nvous pouvez utiliser $lsor ;  par exemple avec Rbase :\n" ,;
  print "  source(\"statgh.r\")\n" ,;
  print "  tdv <- read.table(\"lois.sor\",header=F)\n";
  print "  histProba(tdv[,2], tdv[,1], \" loi...\")\n";
  print "\n" ;

} ; # fin sub decritLoi

## -- fin de lois.pl