Numeric error propagation: Difference between revisions

Content added Content deleted
m (syntax highlighting fixup automation)
(→‎{{header|Perl}}: eliminated 'uninitialized' warnings, added subroutine signatures, additional tweaks)
Line 1,692: Line 1,692:


=={{header|Perl}}==
=={{header|Perl}}==
Following code keeps track of covariance between variables. Each variable with error contains its mean value and components of error source from a set of indepentent variables. It's more than what the task requires.
Following code keeps track of covariance between variables. Each variable with error contains its mean value and components of error source from a set of independent variables. It's more than what the task requires.
<syntaxhighlight lang="perl">use utf8;
<syntaxhighlight lang="perl">use v5.36;
package ErrVar;
use strict;


package ErrVar;
# helper function, apply f to pairs (a, b) from listX and listY
sub zip(&$$) {
my ($f, $x, $y) = @_;
my $l = $#$x;
if ($l < $#$y) { $l = $#$y };


# helper function, apply function 'f' to pairs (a, b) from listX and listY
my @out;
for (0 .. $l) {
sub zip ($f, $x, $y) {
my @out;
local $a = $x->[$_];
$y = [(0) x @$x] unless @$y; # if not defined
local $b = $y->[$_];
push @out, $f->();
push @out, $f->($x->[$_], $y->[$_]) for 0 .. $#$x;
\@out
}
\@out
}
}


use overload
use overload
'""' => \&_str,
'""' => \&_str,
'+' => \&_add,
'+' => \&_add,
'-' => \&_sub,
'-' => \&_sub,
'*' => \&_mul,
'*' => \&_mul,
'/' => \&_div,
'/' => \&_div,
'bool' => \&_bool,
'bool' => \&_bool,
'<=>' => \&_ncmp,
'<=>' => \&_ncmp,
'neg' => \&_neg,
'neg' => \&_neg,
'sqrt' => \&_sqrt,
'sqrt' => \&_sqrt,
'log' => \&_log,
'log' => \&_log,
'exp' => \&_exp,
'exp' => \&_exp,
'**' => \&_pow,
'**' => \&_pow,
;
;


# make a variable with mean value and a list of coefficient to
# make a variable with mean value and a list of coefficient to
# variables providing independent errors
# variables providing independent errors
sub make {
sub make ($x, @v) { bless [$x, @v] }
my $x = shift;
bless [$x, [@{+shift}]]
}

sub _str { sprintf "%g±%.3g", $_[0][0], sigma($_[0]) }


# mean value of the var, or just the input if it's not of this class
# mean value of the var, or just the input if it's not of this class
sub mean ($x) { ref $x && $x->isa(__PACKAGE__) ? $x->[0] : $x }
sub mean {
my $x = shift;
ref($x) && $x->isa(__PACKAGE__) ? $x->[0] : $x
}


# return variance index array
# return variance index array
sub vlist ($x) { ref $x && $x->isa(__PACKAGE__) ? $x->[1] : [] }
sub vlist {
my $x = shift;
ref($x) && $x->isa(__PACKAGE__) ? $x->[1] : [];
}


sub variance {
sub variance ($x) {
return 0 unless ref($x) and $x->isa(__PACKAGE__);
my $x = shift;
my $s;
return 0 unless ref($x) and $x->isa(__PACKAGE__);
$s += $_ * $_ for @{$x->[1]};
my $s;
$s
$s += $_ * $_ for (@{$x->[1]});
$s
}
}


sub covariance {
sub covariance ($x, $y) {
return 0 unless ref($x) && $x->isa(__PACKAGE__);
my ($x, $y) = @_;
return 0 unless ref($x) && $x->isa(__PACKAGE__);
return 0 unless ref($y) && $y->isa(__PACKAGE__);
my $s;
return 0 unless ref($y) && $y->isa(__PACKAGE__);
zip sub ($a,$b) { $s += $a * $b }, vlist($x), vlist($y);

my $s;
$s
zip { $s += $a * $b } vlist($x), vlist($y);
$s
}
}


sub sigma { sqrt variance(shift) }
sub sigma ($v) { sqrt variance $v }


# to determine if a var is probably zero. we use 1σ here
# to determine if a var is probably zero. we use 1σ here
sub _bool {
sub _bool ($x, $, $) {
abs(mean $x) > sigma $x
my $x = shift;
return abs(mean($x)) > sigma($x);
}
}


sub _ncmp {
sub _ncmp ($a, $b, $) {
my $x = shift() - shift() or return 0;
return 0 unless my $x = $a - $b;
return mean($x) > 0 ? 1 : -1;
mean($x) > 0 ? 1 : -1
}
}


sub _neg {
sub _neg ($x, $, $) {
bless [ -mean($x), [map(-$_, @{vlist $x}) ] ];
my $x = shift;
bless [ -mean($x), [map(-$_, @{vlist($x)}) ] ];
}
}


sub _add {
sub _add ($x, $y, $) {
my ($x, $y) = @_;
my ($x0, $y0) = ( mean($x), mean($y));
my ($x0, $y0) = (mean($x), mean($y));
my ($xv, $yv) = (vlist($x), vlist($y));
bless [$x0 + $y0, zip sub ($a,$b) {$a + $b}, $xv, $yv]
my ($xv, $yv) = (vlist($x), vlist($y));
bless [$x0 + $y0, zip {$a + $b} $xv, $yv];
}
}


sub _sub {
sub _sub ($x, $y, $) {
my ($x, $y, $swap) = @_;
my ($x0, $y0) = ( mean($x), mean($y));
if ($swap) { ($x, $y) = ($y, $x) }
my ($xv, $yv) = (vlist($x), vlist($y));
my ($x0, $y0) = (mean($x), mean($y));
bless [$x0 - $y0, zip sub ($a,$b) {$a - $b}, $xv, $yv]
my ($xv, $yv) = (vlist($x), vlist($y));
bless [$x0 - $y0, zip {$a - $b} $xv, $yv];
}
}


sub _mul {
sub _mul ($x, $y, $) {
my ($x, $y) = @_;
my ($x0, $y0) = ( mean($x), mean($y));
my ($x0, $y0) = (mean($x), mean($y));
my ($xv, $yv) = (vlist($x), vlist($y));
my ($xv, $yv) = (vlist($x), vlist($y));
$xv = [ map($y0 * $_, @$xv) ];
$yv = [ map($x0 * $_, @$yv) ];

$xv = [ map($y0 * $_, @$xv) ];
bless [$x0 * $y0, zip sub ($a,$b) {$a + $b}, $xv, $yv]
$yv = [ map($x0 * $_, @$yv) ];

bless [$x0 * $y0, zip {$a + $b} $xv, $yv];
}
}


sub _div {
sub _div ($x, $y, $) {
my ($x, $y, $swap) = @_;
my ($x0, $y0) = ( mean($x), mean($y));
if ($swap) { ($x, $y) = ($y, $x) }
my ($xv, $yv) = (vlist($x), vlist($y));
$xv = [ map($_/$y0, @$xv) ];

my ($x0, $y0) = (mean($x), mean($y));
$yv = [ map($x0 * $_/$y0/$y0, @$yv) ];
bless [$x0 / $y0, zip sub ($a,$b) {$a + $b}, $xv, $yv]
my ($xv, $yv) = (vlist($x), vlist($y));

$xv = [ map($_/$y0, @$xv) ];
$yv = [ map($x0 * $_/$y0/$y0, @$yv) ];

bless [$x0 / $y0, zip {$a + $b} $xv, $yv];
}
}


sub _sqrt {
sub _sqrt ($x, $, $) {
my ($x0, $xv) = ( mean($x), vlist($x) );
my $x = shift;
my $x0 = mean($x);
$x0 = sqrt($x0);
my $xv = vlist($x);
$xv = [ map($_ / 2 / $x0, @$xv) ];
bless [$x0, $xv]
$x0 = sqrt($x0);
$xv = [ map($_ / 2 / $x0, @$xv) ];
bless [$x0, $xv]
}
}


sub _pow {
sub _pow ($x, $y, $) {
my ($x, $y, $swap) = @_;
if ($x < 0) {
die "Can't take pow of negative number $x" if int($y) != $y or $y & 1;
if ($swap) { ($x, $y) = ($y, $x) }
if ($x < 0) {
$x = -$x;
}
if (int($y) != $y || ($y & 1)) {
exp($y * log $x)
die "Can't take pow of negative number $x";
}
$x = -$x;
}
exp($y * log $x)
}
}


sub _exp {
sub _exp ($x, $, $) {
my ($x0, $xv) = ( exp(mean($x)), vlist($x) );
my $x = shift;
bless [ $x0, [map($x0 * $_, @$xv) ] ]
my $x0 = exp(mean($x));
my $xv = vlist($x);
bless [ $x0, [map($x0 * $_, @$xv) ] ]
}
}


sub _log {
sub _log ($x, $, $) {
my ($x0, $xv) = ( mean($x), vlist($x) );
my $x = shift;
bless [ log($x0), [ map($_ / $x0, @$xv) ] ]
my $x0 = mean($x);
my $xv = vlist($x);
bless [ log($x0), [ map($_ / $x0, @$xv) ] ]
}
}


sub _str { sprintf '%g±%.3g', $_[0][0], sigma($_[0]) }
"If this package were to be in its own file, you need some truth value to end it like this.";


package main;
package main;
Line 1,873: Line 1,828:


my $z1 = sqrt(($x1 - $x2) ** 2 + ($y1 - $y2) ** 2);
my $z1 = sqrt(($x1 - $x2) ** 2 + ($y1 - $y2) ** 2);
print "distance: $z1\n\n";
say "distance: $z1";


# this is not for task requirement
# this is not for task requirement
my $a = $x1 + $x2;
my $a = $x1 + $x2;
my $b = $y1 - 2 * $x2;
my $b = $y1 - 2 * $x2;
print "covariance between $a and $b: ", $a->covariance($b), "\n";</syntaxhighlight>output<syntaxhighlight lang="text">distance: 111.803±2.49
say "covariance between $a and $b: ", $a->covariance($b);</syntaxhighlight>
{{out}}

<pre>distance: 111.803±2.49
covariance between 300±2.46 and -350±4.56: -9.68</syntaxhighlight>
covariance between 300±2.46 and -350±4.56: -9.68</pre>


=={{header|Phix}}==
=={{header|Phix}}==
Retrieved from "https://rosettacode.org/wiki/Numeric_error_propagation"