Simple Perl script to generate and fit noisy linear data

fit-data.pl

use strict;
use warnings;

=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

my $filename = $ARGV[0]
	or die "You must tell me where to retrieve the data\n";

# Open the file
open my $in_fh, '<', $filename;

# Pull in the slope and intercept:
my $slope = <$in_fh>;
$slope =~ s/.*: (.*)\n/$1/;
my $intercept = <$in_fh>;
$intercept =~ s/.*: (.*)\n/$1/;

# Pull in the rest of the data:
my @data = <$in_fh>;

# Calculate the sums:
my ($S_x, $S_y, $S_xy, $S_xx);
for (my $x = 0; $x < @data; $x++) {
	my $y = $data[$x];
	$S_x += $x;
	$S_y += $y;
	$S_xy += $x * $y;
	$S_xx += $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";
}