Cpプロットを書くプログラムを出力させるPerlスクリプト

Cpプロットを書くときに、説明変数の数がおおくなってくると、計算しなければならないCp値の数がとんでもないことになってしまいます。説明変数の数が6つくらいでも57個くらいか。まだまだ、増えると手に負いきれなくなるのでプログラムを書いてどうにかしましょう。僕は下のようなPerlのプログラムを書くことで楽をすることにしました。ま、データ解析の本質とは全く関係ないんだけれども。

use strict;
use warnings;
#CpプロットをRに書かせるためのプログラムをはかせるプログラム
sub combinatorial {
    my $m = shift;## m 個を選択する
	my (@array) = @_;## @_ から
	combinatorial_core($#array+1, $m, @array);
}

sub combinatorial_core {
    my $n = shift;## n 個の要素から
	my $m = shift;## m 個を選択する
	## m 個選択済みなら再帰処理を終了
	return @_ if($m==1);
    my @result;
    ## 要素群 P(x) から要素をどこから取り出すか？ 1 〜 n 個のパターン
    for (my $i = 0; $i < @_; $i++) {
	my @data = @_;
	## 要素群 P(x) から i 個目の要素を１つ取り出す
	my $item = splice(@data, $i, 1);
	## 使用済みの要素は破棄する
	@data = splice(@data, $i);
	## 残った要素群 P(x+1) から次の要素を取り出す処理を再起処理にまかせる
	push(@result, map { "$item,$_" } combinatorial_core($#data, $m-1, @data));
    }
    ## 要素群 P(x) の全組み合わせを返す
    return @result;
}

#階乗の計算をさせる。再起で書けよと思ったのは内緒。
sub factorial{
    my $f=1;
    my $n=shift;
    for (my $i = 1; $i <= $n; $i++) {
	$f *= $i;
    }
    return $f;
}

#組合せを計算させる。定義どおり。
sub num_of_comb{
    my ($n, $k) = @_;
    return (factorial($n)) / ( (factorial($k)) * (factorial($n-$k)) );
}

#Cp値の大体の部分の出力をさせるところ
sub print_result{
    my $num = shift;
    my @data = @_;
    my @results = combinatorial($num, @data);
    my $count = 1;
    foreach my $result(@results){
	my @keyarray = split(',',$result);
	print "rss.$num$count<-sum(lsfit(cbind(";
	print join(',',@keyarray);
	print '),log.mid.price)$residuals^2)',"\n";

	print "Cp.",scalar(@keyarray),$count,"<-(rss.",scalar(@keyarray),$count,"/unbiased.sigma.square.hat)-(length(log.mid.price)-2*",(scalar(@keyarray)+1),")\n";
	$count++;
    }
}

#フルモデルに組み込む説明変数をここの配列に入れる。
my @setumei=("mhw","mid.displacement","mid.mpg","is.eu","is.hybrid","is.suv");

print 'unbiased.sigma.square.hat<-sum(lsfit(cbind(',join(',',@setumei),'),log.mid.price)$residuals^2)',"/(length(log.mid.price)-",scalar(@setumei)+1,")","\n";

#決定係数とCp値の値を出力するところ
for(my $i=2; $i<=scalar(@setumei); $i++){
    print_result($i,@setumei);
    print "\n";
}

#Cp.criteriaを出力
print "Cp.criteria<-c(";

for(my $i=2; $i < scalar(@setumei)+1; $i++){
    for(my $j=1; $j < num_of_comb(scalar(@setumei),$i)+1; $j++){
	unless($i==scalar(@setumei) && $j == num_of_comb(scalar(@setumei),$i)){
	    print "Cp.",$i,$j,",";
	}
	#最後のところだけは「,」を付ける必要がない
	else{
	    print "Cp.",$i,$j;
	}
    }
}

print ")\n";

print "number.of.variables<-c(";

for(my $i = 2; $i < scalar(@setumei)+1 ;$i++){
    for(my $j = 0; $j < num_of_comb(scalar(@setumei),$i); $j++){
	unless($i == scalar(@setumei) && $j == num_of_comb(scalar(@setumei),$i)-1){
	    print $i,",";
	}
	else{
	    print $i;
	}
    }
}

print ")\n";

print 'plot(number.of.variables,Cp.criteria,xlab="number of varoables, p (including constant)",ylab="Cp values")',"\n";

print "abline(0,1)\n";

で、こいつをPerlスクリプトとして実行させると*1、下のようなCpプロットがいっぱつで書ける。一件落着。

*1:実行のさせかたくらは調べてねw