#!/usr/bin/perl

use 5.038.2;

# sample data of two series of wages to detect differences which could qualify as discriminant

# Wilcoxon signed rank method

# v.2 - see anotations at the end
# VERSION SUBRUTINA MUESTRAL

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

ini:
# simulación de datos

my $n = 60;
my (@aa, @bb);
for my $i (0..$n-1){
    $aa[$i] = 100*rand();
    $bb[$i] = 100*rand() - (rand()-0.5) *20;
}

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

# usign signatured subs    
# i.e. &wilcoxon2( @x , @y, $n );      # $n same on both, error triggered 

wilcoxon2( \@aa, \@bb, $n );


# 

exit 3;


sub wilcoxon2 {
    
    my ($a,$b,$n) = @_;
    my @aa = @$a;
    my @bb = @$b;
    
    # step 1. Calculate the differences
    my @d;
    $d[$_] = $aa[$_]-$bb[$_] for (0..$n-1);

    # step 2. Rank the absolute values of the differences (|Di|) from smallest to largest, assigning an integer rank value (Ri) to each difference
    
    my @sd = sort { abs($a) <=> abs($b) } @d;
    my %h;
    my $c = 0;
    for my $i (@sd){            # @sd preserves the sign
	$c++;
	$h{$i}=$c;
    }
    
    # step 3. Calculate the sum of the ranks for the positive differences (T+) and the negative differences (T-) separately
    my ( $sumpos, $sumneg ) = (0,0); 
    for my $i (@sd){
	if ($i > 0){
	    $sumpos += $h{$i};
	}elsif ($i < 0){
	    $sumneg += $h{$i};
	}			
    }

    # step 4. Use the smaller of T+ and T- as your Wilcoxon signed-rank statistic (W)
    my $W = ($sumpos < $sumneg) ? $sumpos : $sumneg;

    # step 5. Compare the calculated W with the critical value from the Wilcoxon signed-rank test table, based on your chosen significance level (α) and sample size
    say "W=$W\n";
    say "Search the critical value for W in the Wilconson signed rank table.";
    say "For n=$n and alpha= 0.01, 0.05 or 0.10. If W is higher, the ";
    say "hypothesis of significant differences happen.\n";

    for my $alpha (0.10, 0.05, 0.01){
	print "alpha = $alpha\t\t";
	my $cv = int((($n * ($n + 1)) / 2) * (1 - $alpha / 2));
	print "$cv\n";
    }

    my $z = ($W-($n*($n+1)/4)) / sqrt($n*($n+1)*(2*$n+1)/24);
    say "\nAlternatively, search z/2 in a normal table, where z is = ";
    say "z=$z // Remember to use a two tail normal table.";
    if ($z>0){
	say "P[ ",1-$z/2," \< Z \< ",1+$z/2," ] = ", probnormal2tails($z/2);
    }else{
	say "P[ ",1+$z/2," \< Z \< ",1-$z/2," ] = ", probnormal2tails($z/2);
    }
    say "";
}
    

sub probnormal2tails {

    my $z = shift;
    my $pi = 4*atan2(1,1);
    my $step = 1/10_000_000;
    my $r = 1/2 ;
    
    # densidad de la distribución normal estandarizada
    # Z = (1/sqrt(2*pi))*exp(1)^(-0.5*z*z)
    # my $pi = 4*atan2(1,1) ;
    
    my $static = 1/(10_000_000*sqrt(2*$pi)) ;
    
    my $e = exp(1) ;
    
    for ( my $i = $step; $i <= $z; $i += $step ) {
	    print "\r$i";
	    $r += $static * ($e**(-$i*$i/2));
	}
    say "";
    
    $r = sprintf ("%.8f",$r);  # la precisión no puede ser muy buena
    
#    return $r;
    
#    say "P[ Z <= $z ] = $r";
#    my $rleft = 1 - $r ;
#    say "P[ Z > $z ] = $rleft";
#    $r -= 2 * $rleft;

    return $r;

#    say "P[ -$z <= Z <= $z ] = $r";
    
}    

# exit 2;


__END__
  
Es un test no paramétrico. Hay otros test que sí lo son y pueden
  ofrecer resultados mejores potencialmente.

v.2 La duda es cómo aplicar este análisis no parámetrico tan poderoso
  a una señal y su contador o índice de tiempo. Es decir, una señal
  que tiene ruido seguro (mucha información) a una de intervalos y,
  que por necesidad, va a contener muy poca información. Shannon. El
  primer paso sería hacer de todo ello una subrutina muestral a la que
  se le pasan dos listas.
  
Un ejemplo sería la aplicación a 'sens3', muestras de 60 minutos de 
  actividad media de 5 minutos, la señal, pero la duda es entonces 
  que como sólo hay 12 datos en cada grabación horaria este test 
  debería hacerse con 12*24 datos = 288 pares de datos, que ya es una
  muestra significativa. La subrutina debería devolver un informe del
  tipo de la versión 1

Resumiendo, esta versión debe poner énfasis en ser versión subrutina