gnuplot で線形回帰の相関係数 を計算するには stats を使う

gnuplot でデータの線形回帰計算をする場合は stats コマンドが楽です。ピアソンの相関係数も計算してくれます。 (非線形の場合はfitを使うしかありません。)

(コメントで指摘されたように、相関係数と決定係数とを混同していましたので、訂正してあります。)

usingで変換できるので、簡単な変換で線形化できる場合はstats コマンドが使えます。

stats を使う例

stats コマンドは1列、または2列のデータの統計的なサマリを表示するものです。平均値や最大・最小値などを外部プログラムを使わずに得ることができます。

統計の指標としてピアソンの相関係数も表示されます。

例

次のようなデータを考えます:

1 11
1.5 10.5
2 12
3 13
4 14
5 14.5

このデータファイルを’stats-test.dat’としてstatsコマンドを使うと、以下のような結果が表示されます:

gnuplot> stats 'stats-test.dat'

* FILE: 
  Records:           6
  Out of range:      0
  Invalid:           0
  Column headers:    0
  Blank:             0
  Data Blocks:       1

* COLUMNS:
  Mean:               2.7500           12.5000
  Std Dev:            1.4068            1.4720
  Sample StdDev:      1.5411            1.6125
  Skewness:           0.3367            0.0000
  Kurtosis:           1.7109            1.5000
  Avg Dev:            1.2500            1.3333
  Sum:               16.5000           75.0000
  Sum Sq.:           57.2500          950.5000

  Mean Err.:          0.5743            0.6009
  Std Dev Err.:       0.4061            0.4249
  Skewness Err.:      1.0000            1.0000
  Kurtosis Err.:      2.0000            2.0000

  Minimum:            1.0000 [0]       10.5000 [1]
  Maximum:            5.0000 [5]       14.5000 [5]
  Quartile:           1.5000           11.0000
  Median:             2.5000           12.5000
  Quartile:           4.0000           14.0000

  Linear Model:       y = 1.011 x + 9.721
  Slope:              1.011 +- 0.1356
  Intercept:          9.721 +- 0.4189
  Correlation:        r = 0.9658
  Sum xy:             218.2

gnuplot> 

上記でy = 1.001 x + 9.9721とあるのは、データの1列目をx、2列目をyとして回帰計算をした結果です。

傾き”Slope”, 切片”Intercept”, 相関関数”Correlation”などは、それぞれ変数STATS_slope, STATS_intercept, STATS_correlationとして参照できます。

statsコマンドで回帰計算した結果を参照して回帰直線を描画する例を下に示します:

(Excelの回帰で表示されるR2は決定係数ですので、相関係数とは異なります。 また念のため書いておきますが、決定係数は相関係数の2乗で定義されるものではありません。参考ページ)

unset grid
set terminal pngcairo color enhanced font "Helvetica,18"
set output "fit-test.png"

set xlabel 'x' font "Helvetica,18"
set ylabel 'y' font "Helvetica,18"

set size ratio 0.8
set tics font "Helvetica,18"
set xrange [:]
set yrange [:]
set xtics 1
set ytics 1
set pointsize 2

stats 'fit-test.dat'

R = sprintf("correlation r = %.4f", STATS_correlation )

set label R2 at graph 0.3,0.2
a1 = STATS_intercept
b1 = STATS_slope

f1(x) = a1 + b1 * x

plot 'fit-test.dat' using ($1):($2) title "" with points \
     lc rgb hsv2rgb(0.9, 1, 1) pt 6 ps 2,\
     f1(x) title "" with lines lt 1 lw 2 lc rgb hsv2rgb(0.7,1,1)

unset output

stats で using を使う例

stats は using を使うことができます。上のstats.datの2列目だけにstatsコマンドを適用すると以下のような結果になります:

gnuplot> stats 'stats-test.dat' using 2

* FILE: 
  Records:           6
  Out of range:      0
  Invalid:           0
  Column headers:    0
  Blank:             0
  Data Blocks:       1

* COLUMN: 
  Mean:              12.5000
  Std Dev:            1.4720
  Sample StdDev:      1.6125
  Skewness:           0.0000
  Kurtosis:           1.5000
  Avg Dev:            1.3333
  Sum:               75.0000
  Sum Sq.:          950.5000

  Mean Err.:          0.6009
  Std Dev Err.:       0.4249
  Skewness Err.:      1.0000
  Kurtosis Err.:      2.0000

  Minimum:           10.5000 [1]
  Maximum:           14.5000 [5]
  Quartile:          11.0000 
  Median:            12.5000 
  Quartile:          14.0000 

gnuplot> 

using で列のデータの演算もできます。下の例は2列目を二乗した場合です:

gnuplot> stats 'stats-test.dat' using ($2 ** 2)

* FILE: 
  Records:           6
  Out of range:      0
  Invalid:           0
  Column headers:    0
  Blank:             0
  Data Blocks:       1

* COLUMN: 
  Mean:             158.4167
  Std Dev:           36.8309
  Sample StdDev:     40.3462
  Skewness:           0.0881
  Kurtosis:           1.5046
  Avg Dev:           33.3333
  Sum:              950.5000
  Sum Sq.:       158714.1250

  Mean Err.:         15.0361
  Std Dev Err.:      10.6322
  Skewness Err.:      1.0000
  Kurtosis Err.:      2.0000

  Minimum:          110.2500 [1]
  Maximum:          210.2500 [5]
  Quartile:         121.0000 
  Median:           156.5000 
  Quartile:         196.0000 

gnuplot> 

fit コマンドは使えないのか?

fitコマンドの結果として出て来る FIT_WSSR, FIT_NDFなどを用いてピアソンの相関係数計算する例もありますが、調べた感じでは正しい値を得るのは難しいようです。

その原因はFIT_WSSRはFIT_NDFのポイント数しか評価しないからです。本来は全てのデータ点(y)と期待値(ハットy・フィット値)の差を評価する必要があります。