Simple Perl script to generate and fit noisy linear data

fit-data.pl

use strict;
use warnings;

use v5.10;

=pod

To run, say something like this:

 perl fit-data.pl data.dat

where data.dat is the file created using make-data.pl

=cut

die "You must tell me where to retrieve the data\n"
  unless @ARGV;

chomp( my ($slope, $intercept, @data) = <> );
s/.*?:\s*// for ($slope, $intercept);

# Calculate the sums:
my ($S_x, $S_y, $S_xy, $S_xx);
foreach my $y (@data) {
        state $x = 0;
	$S_x += $x;
	$S_y += $y;
	$S_xy += $x * $y;
	$S_xx += $x * $x;
        $x++
}

# Figure out the slope and intercept from the sums:
my $N_data = scalar(@data);
my $b = ($S_x * $S_y - $N_data * $S_xy) / ($S_x**2 - $N_data * $S_xx);
my $a = ($S_y - $b * $S_x) / $N_data;

print "Actual slope: $slope; computed slope: $b\n";
print "Actual intercept: $intercept; computed intercept: $a\n";

make-data.pl

use strict;
use warnings;

=pod

To run, type something like this at the command line:

 perl make-data.pl > data.dat

You can also specify the number of data points to generate:

 perl make-data.pl 50 > data-50.dat

=cut


# Get the number of points from the command line, or default to 40:
my $N_points = $ARGV[0] || 40;

# Choose a random slope between -1 and 1...
my $slope = rand(2) - 1;
# ... and a random intercept between -10 and 10:
my $intercept = rand(20) - 10;
# Make the magnitude of the noise relative to the size of the slope:
my $noise_magnitude = $slope / 2;

# Print out the slope and intercept:
print "# slope: $slope\n";
print "# intercept: $intercept\n";

# Generate some noisy data:
for my $x (1..$N_points) {
	my $noisy_point = $intercept + $slope * $x + 2 * rand($noise_magnitude) - $noise_magnitude;
	print "$noisy_point\n";
}