View difference between Paste ID: ygvEUW3k and XKyfR39V
SHOW: | | - or go back to the newest paste.
1
#!/usr/bin/perl
2
use strict;
3
use warnings;
4
use feature 'say';
5
my $data = [];
6
while (<>) {
7
  chomp;
8-
  my ($x, $dummy1, $y, $dummy2) = split;
8+
  my ($x, $xe, $y, $ye) = split;
9-
  push @$data, [$x, $y];
9+
  push @$data, [$x, $xe, $y];
10
}
11-
my $a;
11+
my ($a, $b, $a_err, $b_err);
12-
my $b;
12+
13
  my ($sum_w, $sum_xy, $sum_x, $sum_y, $sum_x2) = (0., 0., 0., 0, 0.);
14-
  my $sum_xy = 0.;
14+
15-
  my $sum_x = 0.;
15+
    my ($x, $xe, $y) = @$d;
16-
  my $sum_y = 0.;
16+
    my $w = ($xe == 0.) ? 1 : 1 / $xe;
17-
  my $sum_x2 = 0.;
17+
    $sum_w += $w;
18
    $sum_xy += $w * $x * $y;
19-
    my ($x, $y) = @$d;
19+
    $sum_x += $w * $x;
20-
    $sum_xy += $x * $y;
20+
    $sum_y += $w * $y;
21-
    $sum_x += $x;
21+
    $sum_x2 += $w * $x ** 2;
22-
    $sum_y += $y;
22+
23-
    $sum_x2 += $x ** 2;
23+
24
  my $denom = $sum_w * $sum_x2 - $sum_x ** 2;
25
  $a = ($n * $sum_xy - $sum_x * $sum_y) / $denom;
26-
  $a = ($n * $sum_xy - $sum_x * $sum_y) / ($n * $sum_x2 - $sum_x ** 2);
26+
  $b = ($sum_x2 * $sum_y - $sum_xy * $sum_x) / $denom;
27-
  $b = ($sum_x2 * $sum_y - $sum_xy * $sum_x) / ($n * $sum_x2 - $sum_x ** 2);
27+
  $a_err = $sum_x2 / $denom;
28-
  my $ys = 0.;
28+
  $b_err = $sum_w / $denom;
29
  my $diff = 0.;
30-
    my ($x, $y) = @$d;
30+
31-
    $ys += ($y - $a * $x - $b) ** 2;
31+
    my ($x,, $y) = @$d;
32
    $diff += ($y - $a * $x - $b) ** 2;
33-
  my $error = 3 * sqrt($ys / ($n - 1));
33+
34
  my $error = 3 * sqrt($diff / ($n - 1));
35
  my $i = 0;
36-
    my ($x, $y) = @$d;
36+
37
    my ($x,, $y) = @$d;
38
    splice @$data, $i, 1 if abs($y - $a * $x - $b) >= $error;
39
    $i++;
40
  }
41
  last if $n == @$data;
42-
say "y = ${a} * x + ${b}";
42+
43
say "y = ${a} * x + ${b}";
44
say "a_err: ${a_err}, b_err: ${b_err}";