LeastSquaresMethodWithError

#!/usr/bin/perl
use strict;
use warnings;
use feature 'say';
my $data = [];
while (<>) {
  chomp;
  my ($x, $xe, $y, $ye) = split;
  push @$data, [$x, $xe, $y];
}
my ($a, $b, $a_err, $b_err);
while (1) {
  my ($sum_w, $sum_xy, $sum_x, $sum_y, $sum_x2) = (0., 0., 0., 0, 0.);
  foreach my $d (@$data) {
    my ($x, $xe, $y) = @$d;
    my $w = ($xe == 0.) ? 1 : 1 / $xe;
    $sum_w += $w;
    $sum_xy += $w * $x * $y;
    $sum_x += $w * $x;
    $sum_y += $w * $y;
    $sum_x2 += $w * $x ** 2;
  }
  my $n = @$data;
  my $denom = $sum_w * $sum_x2 - $sum_x ** 2;
  $a = ($n * $sum_xy - $sum_x * $sum_y) / $denom;
  $b = ($sum_x2 * $sum_y - $sum_xy * $sum_x) / $denom;
  $a_err = $sum_x2 / $denom;
  $b_err = $sum_w / $denom;
  my $diff = 0.;
  foreach my $d (@$data) {
    my ($x,, $y) = @$d;
    $diff += ($y - $a * $x - $b) ** 2;
  }
  my $error = 3 * sqrt($diff / ($n - 1));
  my $i = 0;
  foreach my $d (@$data) {
    my ($x,, $y) = @$d;
    splice @$data, $i, 1 if abs($y - $a * $x - $b) >= $error;
    $i++;
  }
  last if $n == @$data;
}
say "y = ${a} * x + ${b}";
say "a_err: ${a_err}, b_err: ${b_err}";