in reply to (contest) Help analyze PM reputation statistics
Here's the code. (Note: I'm truncating the data for conciseness.)
Now, let's run the above program and see the output:#!/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
Woohoo! It looks like we have a good fit. Converting our fitted model into a Perl function that estimates the count of nodes with a given XP, we get the following: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
(Because I fitted the model against log10(count), we had to exponentiate the resulting formula to get an estimation function for count.)sub estimate_count_from_xp($) { my $xp = shift; 10 ** ( 4.442 - 4.194e-2 * $xp + 1.467e-4 * $xp**2 ); }
Just to see how good our model is, take a look at this plot comparing the actual values (dots) versus the estimated values (line). That's pretty much "on the money."
Cheers,
Tom
Tom Moertel : Blog / Talks / LectroTest / PXSL / Coffee / Movie Rating Decoder
|
---|
Replies are listed 'Best First'. | |
---|---|
Re: OK, here's your analysis (w/ picture!)
by zby (Vicar) on Sep 15, 2004 at 07:08 UTC | |
by tmoertel (Chaplain) on Sep 15, 2004 at 07:32 UTC | |
by zby (Vicar) on Sep 15, 2004 at 12:50 UTC |