#!/usr/bin/perl -wl

use strict;
use File::Temp qw( tempfile );

my %rep_stats=(1=>24694, 2=>23551, 0=>22855, ... -41=>1, );

# only keep 0 < XP < 100 because we want the more mainstream
# values and not the far-out ones

my @xp_sorted = grep { $_>0 && $_<100 } sort { $a <=> $b }
                keys %rep_stats;

# generate our R commands

(my $r_commands = <<EOF) =~ s/^    //mg;
    xp <- scan()
    @xp_sorted

    count <- scan()
    @{[ @rep_stats{@xp_sorted} ]}

    summary(lm(log10(count) ~ xp + I(xp^2)))
EOF

# now, we stuff our R commands into a tempfile,
# which we'll use as STDIN

my $tmp = tempfile()
    or die "can't open tempfile: $!";
print $tmp $r_commands;
seek $tmp, 0, 0 or die "can't seek to BOF: $!";
open STDIN, ">&", $tmp or die "can't dup tmp->STDIN: $!";

# finally, we exec R, which will read our commands
# from STDIN (the temp file will be deleted automatically
# when the program exits)

my @cmd = qw(R --no-save --no-init-file --no-restore-data --slave);
exec @cmd;
die "couldn't exec @cmd : $1";  # should never get here

##</code><code>##

Read 99 items
Read 99 items

Call:
lm(formula = log10(count) ~ xp + I(xp^2))

Residuals:
      Min        1Q    Median        3Q       Max 
-0.111459 -0.032358  0.004959  0.024387  0.109470 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.442e+00  1.239e-02  358.49   <2e-16 ***
xp          -4.194e-02  5.719e-04  -73.33   <2e-16 ***
I(xp^2)      1.467e-04  5.541e-06   26.48   <2e-16 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 0.04027 on 96 degrees of freedom
Multiple R-Squared: 0.9975,	Adjusted R-squared: 0.9974 
F-statistic: 1.889e+04 on 2 and 96 DF,  p-value: < 2.2e-16 

##</code><code>##

sub estimate_count_from_xp($) {
    my $xp = shift;
    10 ** ( 4.442 - 4.194e-2 * $xp + 1.467e-4 * $xp**2 );
}