Math::GSL for perl による線形回帰(最小二乗近似)

更に先程のエントリのオマケ。

GSL ( GNU Scientific Library )は大量の計算機算法に対応していますが、線形回帰(最小二乗近似)も、その中の一つ。

前準備

線形回帰=最小二乗近似 を実装する上で、分散、共分散、相関係数 等の公式は、身につけておくべきかと思いますので、メモ

分散

以前のエントリに記載しています
分散の2種類の計算方法 - end0tknrのkipple - web写経開発

共分散

相関係数

更に別の以前のエントリに記載していますが、分散&共分散の式?により構成されます。
perlのStatistics::Basicでピアソンの相関係数を算出 - end0tknrのkipple - web写経開発

分散共分散行列(共分散行列)

単に分散と共分散を列挙したものですが、以下の通りです。

ちなみに3x3の共分散行列は次の通り。

Math::GSL for perl による線形回帰(最小二乗近似)

今回は、定数項のあるモデル( Y=c0+C1*X )と、定数項のないモデル( Y=C1*X )を、それぞれGSLのgsl_fit_linear(), gsl_fit_mul() で算出します。

参考url:
Math::GSL::Fit - search.cpan.org
GSL - GNU Scientific Library - GNU Project - Free Software Foundation
GSL reference manual Japanese translation

#!/usr/local/bin/perl
use strict;
use warnings;
use utf8;
use Encode;
use Data::Dumper;

# http://search.cpan.org/perldoc?Math%3A%3AGSL%3A%3AFit
# http://www.gnu.org/software/gsl/
# http://cbrc3.cbrc.jp/~tominaga/translations/gsl/index.html

use Math::GSL::Fit qw/:all/;


main();

sub main {

    my @norris_x =
        qw/0.2    337.4  118.2  884.6  10.1  226.5  666.3  996.3
           448.6  777.0  558.2  0.4    0.6   775.5  666.9  338.0
           447.5  11.6   556.0  228.1  995.8 887.6  120.2  0.3
           0.3    556.8  339.1  887.2  999.0 779.0  11.1   118.3
           229.2  669.1  448.9  0.5/;
    my @norris_y =
        qw/0.1    338.8  118.1  888.0  9.2   228.1  668.5  998.5
           449.1  778.9  559.2  0.3    0.1   778.1  668.8  339.3
           448.9  10.8   557.7  228.3  998.0 888.8  119.6  0.3
           0.6    557.6  339.3  888.0  998.5  778.9  10.2  117.6
           228.9  668.4  449.2  0.2/;

    #### MODEL: Y = c0 + c1*X の線形回帰
    my ($stat,                          #成功,失敗
        $model_params,                  #モデルのパラメータ
        $covariance,                    #分散共分散行列
        $sum_squared_covariance) =      #残差平方和
            linear_regression_c0_c1_x(\@norris_x, \@norris_y);

    print "MODEL PARAMS(Y = c0 + c1*X) : \n";
    print " c0=$model_params->[0] , c1=$model_params->[1]\n";
    print "COVARIANCE(cov00, cov01, cov11):\n";
    print " cov00=$covariance->[0] , cov01=$covariance->[1] , ";
    print "cov11=$covariance->[1]\n";
    print "SUM OF SQUARED RESIDUALS:\n";
    print " $sum_squared_covariance\n\n";


    #### MODEL: Y = c1*X の線形回帰
    my ($stat_1,                          #成功,失敗
        $model_params_1,                  #モデルのパラメータ
        $covariance_1,                    #分散共分散行列
        $sum_squared_covariance_1) =      #残差平方和
            linear_regression_c1_x(\@norris_x, \@norris_y);
    print "MODEL PARAMS(Y = c1*X) : \n";
    print " c1=$model_params_1->[0]\n";
    print "COVARIANCE(cov11):\n";
    print " cov11=$covariance_1->[0]\n";
    print "SUM OF SQUARED RESIDUALS:\n";
    print " $sum_squared_covariance_1\n\n";
}

# MODEL: Y = c0 + c1*X
sub linear_regression_c0_c1_x {
    my ($norris_x, $norris_y) = @_;

    my ($stat,                          #成功,失敗
        $c0, $c1,                       #モデルパラメータ
        $cov00, $cov01, $cov11,         #分散共分散行列
        $sumsq) =                       #残差平方和
            gsl_fit_linear($norris_x,
                           1,               #stride
                           $norris_y,
                           1,               #stride
                           scalar(@$norris_x) ); #data size
    
    return ($stat, [$c0, $c1], [$cov00, $cov01, $cov11], $sumsq);
}


# MODEL: Y = c1*X
sub linear_regression_c1_x {
    my ($norris_x, $norris_y) = @_;

    my ($stat,          #成功,失敗
        $c1,            #モデルパラメータ
        $cov11,         #c1の分散
        $sumsq) =       #残差平方和
            gsl_fit_mul( $norris_x,
                         1,
                         $norris_y,
                         1,
                         scalar(@$norris_x) );
    return ($stat, [$c1], [$cov11], $sumsq);
}

1;
__END__

↑と書くと、↓になります

$ ./foo.pl 
MODEL PARAMS(Y = c0 + c1*X) : 
 c0=-0.262323073773928 , c1=1.00211681802045
COVARIANCE(cov00, cov01, cov11):
 cov00=0.0542043302231086 , cov01=-7.74327536315676e-05 , cov11=-7.74327536315676e-05
SUM OF SQUARED RESIDUALS:
 26.6173985294235

MODEL PARAMS(Y = c1*X) : 
 c1=1.00174208046979
COVARIANCE(cov11):
 cov11=7.46806595658454e-08
SUM OF SQUARED RESIDUALS:
 27.6112596299331