#!/usr/bin/perl ######################################## #Author: Jesse Lovegren 29-May-2011 #This script generates Rom critical values #for a given alpha and number of observations #based on equations given in Biometrika 77(3):663-5 # # # Usage: rom.pl [n] [alpha] # ######################################## sub fac{ my \$out = 1; my \$i = shift @_; for \$j (1 .. \$i){ \$out *= \$j; } return \$out; } sub combine{ my \$n = shift @_; my \$k = shift @_; my \$out; \$out = &fac(\$n) / (&fac(\$n-\$k) * &fac(\$k)); return \$out; } my \$n = shift @ARGV; my \$alpha = shift @ARGV; unless (\$n, \$alpha){ die ("Usage: rom.pl [n] [alpha]\n");} ############################## print "="x50, "\nRom critical values, n=\$n, alpha=\$alpha\n" , "="x50, "\n"; ########################### my @powers; \$powers[1] = \$alpha; for \$i ( 2 .. \$n ){ \$powers[\$i] = \$powers[\$i-1] + \$alpha ** \$i; } my \$xterm = \$powers[\$n-1]; ########################### my @rom; my \$yterm = 0;; ########################### ########################### MAINLOOP : for \$p ( 1 .. \$n ) { if (\$p == 1) { \$rom[\$p] = \$alpha; next; } elsif (\$p == 2) { \$rom[\$p] = \$alpha/2; next; } else { #the calculable cases YLOOP : for \$q ( 3 .. \$p ){ \$yterm += &combine(\$n,\$q-2) * (\$rom[\$q-1]) ** (\$n-\$q+2); } \$rom[\$p] = (\$powers[\$p-1] - \$yterm) / \$p; } } ########################## PRINTLOOP : for \$r ( 1 .. \$#rom ){ print "a(\$r,\$alpha) = \$rom[\$r]\n"; } print "=" x 50 , "\n";