Numeric error propagation: Difference between revisions

:* If   f = a ± c,   or   f = c ± a   then   '''σ_f = σ_a'''

:* If   f = a ± b   then   '''σ_f² = σ_a² + σ_b²'''

:* If   f = ca   or   f = ac       then   '''σ_f = |cσ_a|'''

:* If   f = ab   or   f = a / b   then   '''σ_f² = f²( (σ_a / a)² + (σ_b / b)²)'''

:* If   f = a^c   then   '''σ_f = |fc(σ_a / a)|'''

# Given coordinates and their errors:<br>x1 = 100 &plusmn; 1.1<br>y1 = 50 &plusmn; 1.2<br>x2 = 200 &plusmn; 2.2<br>y2 = 100 &plusmn; 2.3<br> if point p1 is located at (x1, y1) and p2 is at (x2, y2); calculate the distance between the two points using the classic pythagoreanPythagorean formula: <br> '''<big><big> d = &radic;(<span style="text-decoration:overline"> &nbsp; (x1 - x2)²² &nbsp; + &nbsp; (y1 - y2)² &nbsp; <sup/span>2 </supbig></big>)'''

# Print and display both &nbsp; <big> '''d''' </big> &nbsp; and its error.

<langsyntaxhighlight Adalang="ada">generic

end Approximation;</langsyntaxhighlight>

<langsyntaxhighlight Adalang="ada">with Ada.Text_IO;

end Approximation;</langsyntaxhighlight>

<langsyntaxhighlight Adalang="ada">with Approximation, Ada.Numerics.Elementary_Functions;

end Test_Approximations;</langsyntaxhighlight>

<langsyntaxhighlight lang="algol68"># MODE representing a uncertain number #

print( ( "distance: ", fixed( v OF d, 0, 2 ), " +/- ", fixed( uncertainty OF d, 0, 2 ), newline ) )</langsyntaxhighlight>

<langsyntaxhighlight lang="cpp">#pragma once

};</langsyntaxhighlight>

<langsyntaxhighlight lang="cpp">#include <cstdlib>

}</langsyntaxhighlight>

<langsyntaxhighlight lang="d">import std.stdio, std.math, std.string, std.typecons, std.traits;

}</langsyntaxhighlight>

Following the propagation of error estimates through a computation is a nightmare of convoluted formulae wherein mistakes are easily made. The basic method is to derive the formulae according to the standard rules while carrying forward the calculation by hand with the utmost caution. A computer can perform the calculations, but the real problem is in ensuring that it performs the correct calculations, not some incorrectmisbegotten calculationconfusion... So, rather than attempt to "optimise" the calculation, the objective is to reduce brain strain by producing code whose ''checkability'' is optimised instead, somewhat as follows: <langsyntaxhighlight Fortranlang="fortran"> PROGRAM CALCULATE !A distance, with error propagation.

END !Enough.</langsyntaxhighlight>

All outputs are formatted for the specific values of the test data so as to avoid excessive spacing. Since the given errors are of the order of one, showing more than one decimal digit in the result would be silly. Similarly, standard precision suffices though for more elaborate calculations, double precision could well be preferable for the working out.

Rather than agonise over devising adjoint formulae for the error propagation through some calculation, one can perform the desired calculation via routines that will carry along the error term, as follows:<langsyntaxhighlight Fortranlang="fortran"> MODULE ERRORFLOW !Calculate with an error estimate tagging along.

DATAINTEGER VSP/0/VOUT !StartOutput with nothingfile.

STACKV(VSP - 1) = STACKV(VSP - 1) + STACKV(VSP) !Do the deed.

STACKE(VSP - 1) = SQRT(STACKE(VSP - 1)**2 + STACKE(VSP)**2) !The errors follow.

VSP = VSP - 1 !Two values have become one.

VSP = VSP - 1 !Two values have become one.

STACKV(VSP) = STACKV(VSP)**2 !This will never be negative.

INTEGER P !Though only integer powers for this routine, so no EXP(P*LOGLN(x)).

USE ERRORFLOW !For the details.

Calculate SQRT([(X2 - X1)**2 + (Y2 - Y1)**2)]

VSPCALL = 0 VINIT(6) !Start my arithmetic.

WRITE (6,2) STACKV(1),CPM,STACKE(1) !Ahh, the relief.

END !Enough.</langsyntaxhighlight>

This is closer to the idea of extending the language to supply additional facilities. At the cost of hand-compiling the arithmetic expression into a sequence of pseudo-machine code subroutines, it is apparent that the mind-tangling associated error formulae need no longer be worried over. The various arithmetic subroutines have to be coded correctly with careful attention to the V or E statements in fact being for the V and E terms (cut&paste followed by inadequate adjustment the culprit here: such mistakes are less likely when using a card punch because copying is more troublesome), but this is a straightforward matter of checking. And indeed the VPOW(2) routine has the same effect as VSQUARE. The output is the same. Output:

Although the example source file uses the F90 MODULE protocol, this is merely for convenience. Otherwise, there would be the same collection of routines, all sharing a COMMON work area. In ore general use one would extend the routines to provide more checking and some sort of optional trace facility. Indeed, itIt would be easy enough to prepare an interactive calculator using this scheme so that different calculations (and data) could be more easily experimented with. Inspection might suggest a return to the laboratory in order to measure some factor with greater precision because its error proves to be the largest contributor to the spread in the result.

A latter-day expansion of Fortran makes it possible to define a compound entity such as a value and its associated error, for instance,<langsyntaxhighlight Fortranlang="fortran"> TYPE DATUM

Whereupon, instead of a swarm of separate variables named according to some scheme, you can have a collection of variables with subcomponents named systematically. Further, via a great deal of syntax one can devise functions dealing with those compound types and moreover, prepare procedures that will perform operations such as addition and subtraction, etc. and merge these with the ordinary usages of addition, etc. of ordinary variables.. One would then write the formula to be calculated, and it would all just happen. This has been done for arithmetic with rational numbers, in [[Arithmetic/Rational#Fortran]] for example.

<langsyntaxhighlight lang="go">package main

}</langsyntaxhighlight>

<langsyntaxhighlight lang="haskell">data Error a = Error {value :: a, uncertainty :: a} deriving (Eq, Show)

<langsyntaxhighlight lang="unicon">record num(val,err)

end</langsyntaxhighlight>

<langsyntaxhighlight lang="j">num=: {."1

dist=: +/&.:*:</langsyntaxhighlight>

<langsyntaxhighlight lang="j">x1=: 100 unc 1.1

y2=: 100 unc 2.3</langsyntaxhighlight>

<langsyntaxhighlight lang="j"> num x1

1.1</langsyntaxhighlight>

<langsyntaxhighlight lang="j"> num 100

0</langsyntaxhighlight>

<syntaxhighlight lang="text"> 3 dist 4

5</langsyntaxhighlight>

<langsyntaxhighlight lang="j">add=: +&num unc dist&unc

exp=: ^&num unc |@(^&num * dist&(unc%num))</langsyntaxhighlight>

<langsyntaxhighlight lang="j"> exp&0.5 (x1 sub x2) add&(exp&2) y1 sub y2

111.803 2.48717</langsyntaxhighlight>

<langsyntaxhighlight lang="java">public class Approx {

Approx x2y1 = new Approx(50, 1.2);

Approx y1x2 = new Approx(200, 2.2);

}</langsyntaxhighlight>

<pre>111.80339887498948±2.9383668933610044871670631463423</pre>

<langsyntaxhighlight lang="scala">import java.lang.Math.*

val x2y1 = Approx(50.0, 1.2)

val y1x2 = Approx(200.0, 2.2)

println(((x1 - x2).pow(2.0) + (y1 - y2).pow(2.0)).pow(0.5)) // => 111.80339887498948 ±2.938366893361004

}</langsyntaxhighlight>

<pre>111.80339887498948 ±2.9383668933610044871670631463423</pre>

=={{header|Mathematica}}/{{header|Wolfram Language}}==

<langsyntaxhighlight lang="mathematica">PlusMinus /: a_ ± σa_ + c_?NumericQ := N[(a + c) ± σa];

PlusMinus /: (a_ ± σa_)^c_?NumericQ := N[a^c ± Abs[a^c σa/a]];</langsyntaxhighlight>

<langsyntaxhighlight lang="mathematica">x1 = 100 ± 1.1;

d = Sqrt[(x1 - x2)^2 + (y1 - y2)^2]</langsyntaxhighlight>

<langsyntaxhighlight lang="parigp">add(a,b)=if(type(a)==type(b), a+b, if(type(a)=="t_VEC",a+[b,0],b+[a,0]));

pow(add(pow(sub(x1,x2),2),pow(sub(y1,y2),2)),.5)</langsyntaxhighlight>

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 indepententindependent variables. It's more than what the task requires.

<langsyntaxhighlight lang="perl">use utf8v5.36;

forsub zip (0$f, ..$x, $ly) {

push @out, $f->($x->[$_], $y->[$_]) for 0 .. $#$x;

'""' => \&_str,

'+' => \&_add,

'-' => \&_sub,

'*' => \&_mul,

'/' => \&_div,

'bool' => \&_bool,

'<=>' => \&_ncmp,

'neg' => \&_neg,

'sqrtlog' => \&_sqrt_log,

'logexp' => \&_log_exp,

'exp**' => \&_exp_pow,

sub make ($x, @v) { bless [$x, @v] }

sub variance ($x) {

sub covariance ($x, $y) {

return 0 unless ref($xy) && $xy->isa(__PACKAGE__);

my $s;

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

sub _bool ($x, $, $) {

sub _ncmp ($a, $b, $) {

return 0 unless my $x = shift()$a - shift() or return 0$b;

return mean($x) > 0 ? 1 : -1;

sub _neg ($x, $, $) {

sub _add ($x, $y, $) {

my ($xx0, $yy0) = @_( mean($x), mean($y));

my ($x0xv, $y0yv) = (meanvlist($x), meanvlist($y));

sub _sub ($x, $y, $) {

my ($xx0, $y,y0) = ( mean($swapx), = @_mean($y));

if ($swap) { my ($xxv, $yyv) = (vlist($yx), vlist($xy) });

my ( bless [$x0, - $y0), =zip sub (mean($xa,$b) {$a - $b}, mean($y));xv, $yv]

sub _mul ($x, $y, $) {

my ($xx0, $yy0) = @_( mean($x), mean($y));

my ($x0xv, $y0yv) = (meanvlist($x), meanvlist($y));

my ( $xv, $yv) = (vlist[ map($x)y0 * $_, vlist(@$y)xv) ];

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

sub _div ($x, $y, $) {

my ($xx0, $y,y0) = ( mean($swapx), = @_mean($y));

if ($swap) { my ($xxv, $yyv) = (vlist($yx), vlist($xy) });

my ($x0, $y0)yv = (mean[ map($x)x0 * $_/$y0/$y0, mean(@$y)yv) ];

sub _sqrt ($x, $, $) {

my $x0 = meansqrt($xx0);

my $xv = vlist[ map($x_ / 2 / $x0, @$xv) ];

sub _pow ($x, $y, $) {

my if ($x, $y,< $swap0) = @_;{

if ( $x <= 0) {-$x;

sub _exp ($x, $, $) {

sub _log ($x, $, $) {

printsay "distance: $z1\n\n";

printsay "covariance between $a and $b: ", $a->covariance($b), "\n";</langsyntaxhighlight>output<lang>distance: 111.803±2.49

<langsyntaxhighlight PicoLisplang="picolisp">(scl 12)

(round (sqrt (cdr N) 1.0) 8) ) ) )</langsyntaxhighlight>

<langsyntaxhighlight PicoLisplang="picolisp">(de distance (X1 Y1 X2 Y2)

(unc 100. 2.3) ) ) )</langsyntaxhighlight>

<langsyntaxhighlight lang="python">from collections import namedtuple

p1, p2, distance(p1, p2)))</langsyntaxhighlight>

<langsyntaxhighlight lang="racket">#lang racket

(norm (±- x1 x2) (±- y1 y2))</langsyntaxhighlight>

<langsyntaxhighlight lang="ruby">class NumberWithUncertainty

puts ((x1 - x2) ** 2 + (y1 - y2) ** 2).sqrt</langsyntaxhighlight>

{{out]]}}

<langsyntaxhighlight lang="scala">import java.lang.Math._

}</langsyntaxhighlight>

<langsyntaxhighlight lang="tcl">package require Tcl 8.6

}</langsyntaxhighlight>

<langsyntaxhighlight lang="tcl">RAII-class create Err {

}</langsyntaxhighlight>

<langsyntaxhighlight lang="tcl">set x1 [Err 100 1.1]

puts "d = [$d p]"</langsyntaxhighlight>