#!/usr/bin/perl -w

use 5.030;

# Logaritmic regression
# y = a + b * log(x)            # (el LOG en BASIC es neperiano)

my (@x, @y, @s);
@s = (0) x 6;
my $n = 0;
while (<DATA>){
    chomp;
    last if (/__END__/);
    $n++;
    ($x[$n], $y[$n]) = split /\,/, $_;
    $x[$n] = log($x[$n]);
    $s[1] += $x[$n];
    $s[2] += $y[$n];
    $s[3] += $x[$n]*$x[$n];
    $s[4] += $y[$n]*$y[$n];
    $s[5] += $x[$n]*$y[$n];
#    $n++;
}
# $n--;



$x[0]=1;
$y[0]=1;
my @comprobacion = logregresion(\@x, \@y);
print "alfa beta R2 => @comprobacion \n";

my $c = $n*$s[5]-$s[2]*$s[1];
my $d = $n*$s[3]-$s[1]*$s[1];
my $b = $c/$d;
my $a = ($s[2]-$b*$s[1])/$n;

say "Modelo de linea ajustada  y = $a + $b * log( x )";
$s[4] += -$s[2]*$s[2] / $n;
$s[1] = $b *( $s[5]-$s[1]*$s[2]/$n );
$s[2] = $s[4]-$s[1];
$s[5] = $s[1] / $s[4];
say "Coeficiente de correlación = ", sqrt($s[5]), "";
say "Coeficiente de determinación = $s[5]";
say "Varianza de la estimación = ", $s[2]/($n-2), "";
say "Desviación típica = ", sqrt( $s[2]/($n-2) ), "";
inq:
say "";
say "Interpolación:";
say "pulse 1: Estimación de y conocida x";
say "pulse 2: Estimación de x conocida y";
say "pulse otra cosa para salir...";

chomp (my $r = <STDIN>);

if ($r eq "1"){
    print "x = ";
    chomp ($r = <STDIN>);
    print "y = ", $a + $b * log($r), "\n";
    goto inq;
}elsif ($r eq "2"){
    print "y = ";
    chomp ($r = <STDIN>);
    print "x = ", exp( ($r-$a)/$b ), "\n";
    goto inq;
}else{
    exit 2;
}

exit 1;

sub log10 {
    return (log(shift)/log(10));
}

sub logregresion {
    my ($u, $v) = @_;
    my @x = @$u;
    my @y = @$v;
    my ($sumxy, $sumx, $sumy, $sumx2);
    for my $i (0 .. scalar(@x)-1) {                                   # 0j0, en la versión log se toma de 1 a $n, no desde 0
	next if $i == 0;
#	$x[$i] = log($x[$i]);     # el log es de x, pero YA ESTA CALCULADO ARRIBA
	$sumxy+=$x[$i]*$y[$i];
	$sumx+= $x[$i];
	$sumy+= $y[$i];
	$sumx2 += $x[$i]*$x[$i];
    }
    my $count = $n;
    my $beta = ($count*$sumx2-$sumx**2 != 0) ? (($count*$sumxy-$sumx*$sumy) / ($count*$sumx2-$sumx**2)) : 1;
    my $alfa = $sumy/$count - $beta*$sumx/$count;
        # print "-" x 70,"\n";
        # print "y = $alfa + $beta x\n";
        # print "-" x 70,"\n";
    my $media = $sumy/$count;
        # print "\tD\tF\tAD\tSE\n";
    
    my ($ra, $rb, $j);
#    $j = -1;                  # para datos x,y de 1 a n
    for my $i (0 .. scalar(@x)-1){
	next if $i == 0;
#	$j++;
#	next if $j == 0;
	#    $dd=$y[$i];
	# $ad=abs($dd-$alfa-$beta*$x[$i]);
	# $se=($dd-$alfa-$beta*$x[$i])**2;
	#        printf "%i)\t%.02f\t%.2f\t%.2f\t%.2f\n",$i,$dd,$f,$ad,$se;
	$ra += ($alfa+$beta*$x[$i]-$media)**2;
	$rb += ($y[$i]-$media)**2;      # originalmente $j-1
	#    $mad+=$ad;
	# $mse+=$se;
    }
    my $R2 = ($rb != 0) ? ($ra/$rb) : 1;
        # $mad /= $count;
        # $mse /= $count;
        # print "R2=$r2\tMAD=$mad\tMSE=$mse\n";
        # for $i ($count+1..$count+10){
        # print "$i)\t\t", $alfa+$beta*$i, "\n";
        # }
    return ($alfa, $beta, $R2);
}




__DATA__
15, 12.1
25, 13.4
35, 14.3
75, 16.3
95, 16.8
125, 17.6
145, 17.9
215, 18.9
__END__

Este código fue tomado del BASIC de J.R. VIZMANOS, 1984.