確率・統計

(11) 二標本の解析 -1-

今まで紹介した推定や検定は、一標本の分布に対する手法が主体となっていました。二標本の分布が同じものであるかを判別する検定法も一部ありましたが、これらはそれぞれの分布の標本間に何の関連性もないことが前提となっています。この章では、二つのデータにある関連性の調べ方や、関連性のあるデータを解析するための手法について紹介したいと思います。

1) 散布図(Scatter Plot)

「相関係数(Correlation Coefficient)」や「回帰曲線(Regression Curve)」の定義については、「(4) 多変数の確率分布」の中で前に紹介していますが、もう少し具体的なデータを使ってもう一度その内容を確認しておきたいと思います。

上に示した表は、統計局ホームページから取得した、人口(千人単位)と年間の交通事故発生件数を各都道府県別に表したものです。集計時の年代に差がありますが、同じ年代のデータが揃わなかったため近い年代のものを使っています。下図は、それぞれの分布をヒストグラムで表したものです。

ある程度予想できる結論にはなりますが、人口の多い都道府県ほど交通事故の発生件数は多くなる傾向がヒストグラムからもはっきりとわかります。発生件数だけを見ると東京が突出しているように見えますが、「人が多いほど事故が多い」と考えれば当然の結論になるわけです。そこで、交通事故の発生件数を人口(千人単位)で割った比率を求めて分布を調べると、下図のような結果になります。

発生件数を千人単位の人口で割ったということは、千人あたりで換算した発生件数と考えることもできるので、10 を超えている群馬・静岡・香川・佐賀などは百人あたり年間一回程度は事故に遭っていることになります。

次に、横軸を人口、縦軸を交通事故発生件数としてグラフに点をプロットしてみます。

人口が多くなるほど交通事故の発生件数が多くなる様子が、このグラフからはっきりと認識できます。このようなグラフは「散布図(Scatter Plot ; Scattergraph)」と呼ばれ、二変数間の関係を確認するときに便利です。複数の変数に対する関連性のことを「相関(Correlation)」といい、散布図を見ると、交通事故発生件数は人口に対して直線状に増加する傾向が見られるので、このような場合は「正の相関がある」ことになります。逆に、横軸方向に増加するに従って縦軸が減少するような場合、「負の相関がある」ことになり、そのどちらとも判別できないようなときは「相関がない」と判断されます。

人口の少ない地域は過疎化のため高齢者の占める割合が多いことが予想されます。そこで、65 歳以上の人口が全人口に対して占める割合を都道府県別に計算すると、次のような結果になります。

都道府県 65歳以上の 人口比率(H.20) 都道府県 65歳以上の 人口比率(H.20) 都道府県 65歳以上の 人口比率(H.20) 北海道 23.58 石川 22.86 岡山 24.28 青森 24.43 福井 24.26 広島 23.00 岩手 26.26 山梨 23.65 山口 26.93 宮城 21.54 長野 25.52 徳島 26.07 秋田 28.43 岐阜 22.90 香川 24.83 山形 26.68 静岡 22.63 愛媛 25.62 福島 24.17 愛知 19.17 高知 27.81 茨城 21.32 三重 23.15 福岡 21.45 栃木 21.08 滋賀 19.69 佐賀 23.95 群馬 22.47 京都 22.37 長崎 25.21 埼玉 19.13 大阪 21.21 熊本 25.10 千葉 20.14 兵庫 22.07 大分 25.92 東京 20.24 奈良 22.65 宮崎 25.18 神奈川 19.23 和歌山 26.09 鹿児島 26.03 新潟 25.55 鳥取 25.55 沖縄 17.22 富山 25.25 島根 28.55

予想通り、人口が少ないほど高齢者の占める割合は増加する傾向にあるようです。このような場合は「負の相関がある」ということになります。

2) 相関係数と回帰曲線

相関図を見れば、二変数間の関連性を読み取ることが可能になります。関連性があると判断できたとき、それがどの程度の強さなのかを定量的に測るための指標として「相関係数(Correlation Coefficient)」があります。多変数の確率密度関数から求められる理論上の相関係数は「母相関係数(Population Correlation Coefficient)」で「相関係数(Correlation Coefficient)」の章で、また、実際のデータから求められる相関係数は「標本相関係数(Sample Correlation Coefficient)」で「主成分分析(Principal Component Analysis ; PCA)」の章でそれぞれ紹介しています。

母相関係数を ρ、標本相関係数を r で表すと、それぞれの式は以下のようになります。

σ_xy は「母共分散(Population Covariance)」 E[( x - μ_x )( y - μ_y )] = ∫_x∈R∫_y∈R( x - μ_x )( y - μ_y ) p( x, y ) dxdy を表します。ここで p( x, y ) は x, y を確率変数とする二変数の確率密度関数、μ_x と μ_y はそれぞれ x, y の母平均であり、σ_xy は ( x - μ_x )( y - μ_y ) に対する p( x, y ) の期待値を意味します。これに対し、s_xy は「標本共分散(Sample Covariance)」であり、確率密度確率 p( x, y ) に従う標本 ( x_i, y_i ) ( i = 1, 2, ... N ) を使った ( x_i - m_x )( y_i - m_y ) の和で求められます。ここで、m_x と m_y はそれぞれ x_i, y_i から得られる標本平均を表しています。
σ_x, σ_y は p( x, y ) に対する x, y の母標準偏差であり、s_x, s_y は標本標準偏差を表します。ちなみに s_xy は N で、s_x, s_y は √N で和を割っているため、分母と分子で N が消されることで r は上記に示した最後の式のようになります。

二つの変数の関連性を定量的に調べるための他の方法として、曲線への当てはめがあります。一般的には、曲線を M 個の関数 f_i(x) ( i = 1, 2, ... M ) による線形結合

で表すことができると仮定して、まずは N 個のデータ ( x_j, y_j ) ( j = 1, 2, ... N ) の x_j を上記関数に代入し、y を求めます。求めた y は、実測値である y_j とは完全に一致することはなく、ある程度のズレがあるはずです。そこで、求めた y と実測値 y_j の差の二乗を各データに対して求め、その和を求めた上で、和が最小となるような係数 a₀, a₁, a₂, ... a_M を求めます。この手法を「最小二乗法(Least Squares)」といいます。

直線へ当てはめる場合は y = ax + b として計算すればよいので、求める和を J とすると、

で求めることができます。2 で割っているのは計算時に生じる係数を打ち消すためで、結果には影響しません。J が最小となるときの係数 a, b は、∂J / ∂a = 0、∂J / ∂b = 0 を解くことで得ることができて、

となります。これは連立一次方程式の形であり、

と表されます。これを「正規方程式」といいます。上側に Σ_j{1→N}( 1 ) = N、下側に Σ_j{1→N}( x_j ) を掛けて辺々引くと

となって、a を求めることができますが、上式の中の a に対する係数は ( データの二乗和 x N ) と ( 和の二乗 ) の差となっていて、データ x の標本平均を m_x としたとき

となることから分かるように、この値は標本分散 v_x = Σ_j{1→N}( ( x_j - m_x )² ) に N を掛けた N・v_x と等しくなります。右辺はデータ y の標本平均を m_y として

より、標本共分散 s_xy = Σ_j{1→N}( ( x_j - m_x )( y_j - m_y ) ) に N を掛けた N・s_xy に等しくなるので、係数 a の値は

と非常にシンプルな式となります。また、この値を二式のうちのどちらかに代入すれば、b の値は

と求められます。よって、直線の方程式は

になります。x = ay + b とした場合も同様の計算によって、

という結果が得られます。どちらの直線も点 ( m_x, m_y ) を通りますが、その傾きは異なります。前者の式は、x に対して y がどのように変化しているかを表し、逆に後者の式は、y に対して x がどのように変化しているかを表しています。そのため、y = ax + b としたときの x は「独立変数(Independent Variable)」または「説明変数(Explanatory Variable)」と呼ばれ、y は「従属変数(Dependent Variable)」または「目的変数(Objective Variable)」と呼ばれます。英語では一般的に "Independent Variable", "Dependent Variable" という名が用いられ、説明変数や目的変数の呼び名はあまり使われないようです。「独立変数」という名前で表されるように、x はそれぞれ独立したデータで、y は各 x に従う値であると考えるとこの名前の意味が理解できるのではないかと思います。
誤差が全く無く、全てのデータが直線上にある状態ならば、y = a( x - m_x ) + m_y に対して x = (1/a)( y - m_y ) + m_x となって、どちらの直線も一致します。これが一致しないのは誤差の見方によるもので、y = a( x - m_x ) + m_y では各データを表す点から y 軸に沿って、当てはめた直線に伸ばした線分の距離が全体的に最も小さくなるようにしているのに対し、x = a( y - m_y ) + m_x の場合は x 軸に沿った線分の距離を最小にする式を求めているという違いから生じます。

前節で紹介した「人口」と「交通事故発生件数」の相関では、人口が多いから交通事故も多くなると考えることができます。よって、人口が独立変数で、交通事故発生件数が従属変数であると考えるのが自然です。人口に比例して交通事故も多くなると考えたとき、誤差によってある程度のズレは発生するものの、それは直線で近似することができます。その式は、人口を x、交通事故発生件数を y としたとき、y = ( s_xy / v_x )( x - m_x ) + m_y で求められることになります。実際に求めてみると、

より y = 5.95x + 1.56 x 10³ になります。

「回帰曲線(Regression Curve)」の定義は、次のような式で表されるのでした。

p(y|x) = p( x, y ) / p_x( x ) は、二変数の確率密度関数 p( x, y ) における y の x に関する条件付確率密度であり、x をある値に固定したときの y の確率密度関数を表しています。よって、μ_y(x) は、それぞれの x について y の平均値を求めてそれをプロットしたときの曲線を表すことになります。

最小二乗法は、回帰曲線を推定するための一つの手法と考えることができます。各 x について y はある確率密度分布に従って分布していると仮定して、点全体の分布が確率的に最も高くなるような曲線を回帰曲線と推定するわけです。このとき、後て述べるように「最尤法(Maximum Likelihood Estimation ; MLE)」が利用されています。

二変数の確率密度関数 p( x, y ) において、回帰曲線 y = μ_y(x) が M 個の関数 f_i(x) ( i = 1, 2, ... M ) による線形結合

で表されると仮定します(ここで、f₀(x) = 1 とします)。N 個の測定値 ( x_j, y_j ) ( j = 1, 2, ... N ) は、真の値 μ_y(x_j) に対してランダムな誤差 ε_j が加わって y_j になると考えれば

と表すことができます。誤差 ε_j が平均 0、分散 σ² の正規分布 N( 0, σ² ) に従い、しかも分散は x によらず常に一定であると仮定すると、尤度 L( ε₁, ε₂, ... ε_N ; a₀, a₁, a₂, ... a_M ) = L( ε, a ) は

となり、対数を取って

になります。尤度を最大にするためには、上式の右辺が最大になればよく、変数となる係数 a_i を持つ項は -Σ_j{1→N}( ε_j ) / 2σ² なので、

が最小になるような a_i を求めればよいことになります。これは、最小二乗法そのものであり、求められた係数 a_i は最尤推定量であるということになります。

直線への当てはめを行うとき、y ( または x ) だけの誤差ではなく、両方の誤差が混在していると考えると、真の値を ( μ_xj, μ_yj ) としたとき、各データ ( x_j, y_j ) は

と表すことができます。ここで、( ε_xj, ε_yj ) は各データに対する x, y 方向の誤差成分を表します。この誤差成分が平均 0、分散 σ² の正規分布 N( 0, σ² ) に従い、しかも分散は x, y によらず常に一定であると仮定すると、尤度 L( ε_x1, ε_x2, ... ε_xN ; ε_y1, ε_y2, ... ε_yN ) = L( ε_x, ε_y ) は

であり、対数を取って

となるので、

が最小になるような ε_xj, ε_yj を求めればよいことになります。

を J の式に代入すれば、

となるので、真の値 ( μ_xj, μ_yj ) が直線 ax + by + c = 0 上にあるとして、右辺が最小になるような ( μ_xj, μ_yj ) を、aμ_xj + bμ_yj + c = 0 という制約条件下で求めると(補足1)、

より a / b = α, c / b = β として J の式に代入すれば

これを μ_xj で偏微分すると

なので、左辺がゼロになる時

になります。μ_yj に対しても同様の手法で求めると

となって、これを J の式に戻すと

ここで、ax + by + c = 0 の両辺に任意のゼロでない数を掛けても直線は変化しないことから、a² + b² = 1 と正規化すれば、

となります。これを c で偏微分すると

なので、∂J / ∂c = 0 のとき

と求められます。ここで、m_x, m_y はそれぞれ x_j, y_j から得られる標本平均を表しています。これをもう一度 J の式に代入すれば

と変形できますが、二次形式

は、対象行列

とベクトル n = ( a, b )^T を使って

と表すことができるのを利用すれば、

としたとき J = ( n, Vn ) / 2 であり、a² + b² = 1 より ||n|| = 1 すなわち n は単位ベクトルということになります。

J = Σ_j{1→N}( ( x_j - μ_xj )² + ( y_j - μ_yj )² ) / 2 という式は、( x_j, y_j ) から ( μ_xj, μ_yj ) までの距離の二乗の和を表しています。( μ_xj, μ_yj ) が直線 ax + by + c = 0 上にある場合、J が最小になるためには 線分 ( x_j, y_j ) - ( μ_xj, μ_yj ) が直線 ax + by + c = 0 に垂直である必要があります。このとき、a² + b² = 1 と正規化すれば J = Σ_j{1→N}( ( ax_j + by_j + c )² ) / 2 という式で表されることになりますが、これは直線 ax + by + c = 0 の式に ( x, y ) = ( x_j, y_j ) を代入した値の二乗和を意味します。さらに J が最小になるように c を定めると、直線は ( x_j, y_j ) の標本平均 ( m_x, m_y ) を通ることになり、この時 J = Σ_j{1→N}( { a( x_j - m_x ) + b( y_j - m_y ) }² ) / 2 になります。

上に示したように、この式は二次形式であり、( n, Vn ) / 2 の形で表すことができます。

とすれば V は

で表されます。J を N で割っても最小値に対する n に変化はないので、V の要素を N で割っても問題はありません。すると、V は「(標本)共分散行列 ( (Sample) Covariance Matrix )」 そのものになります。よって、共分散 E[( x - μ_x )( y - μ_y )] を要素とする「(母)共分散行列」を抽出した標本から推定するとき、誤差が正規分布 N( 0, σ² ) に従うと仮定すれば、最尤推定量として標本共分散行列の V を利用することができて、共分散の最尤推定量は s_xy になることが以上のことから証明されたことになります。

「人口」と「交通事故発生件数」において、v_x = 6.81 x 10⁶, v_y = 2.63 x 10⁸, s_xy = 4.05 x 10⁷ より V は

になります。V の固有値 λ は | λE - V | = ( λ - v_x )( λ - v_y ) - s_xy² の解であり、

で求められ、その解は 2.70 x 10⁸, 5.74 x 10⁵ になります。λ = 2.70 x 10⁸ のときの固有ベクトルは、連立方程式

を解くことで求められ、両式とも b = 6.49a という式に変形できるので、a² + b² = 1 という制約から ( a, b ) = ( 0.152, 0.988 ) という結果が得られます。λ = 5.74 x 10⁵ の場合も同様の方法によって、( a, b ) = ( -0.988, 0.152 ) と求めることができます。この二つの固有ベクトルは互いに直交し、最大固有値である λ = 2.70 x 10⁸ に対する単位ベクトル ( a, b ) = ( 0.152, 0.988 ) は、二次形式 f = ( n, Vn ) で表される楕円体の長軸方向を、最小固有値である λ = 5.74 x 10⁵ に対する単位ベクトル ( a, b ) = ( -0.988, 0.152 ) は楕円体の短径をそれぞれ表すことになります。よって、直線 ax + by + c = 0 は、点 ( m_x, m_y ) = ( 2.72 x 10³, 1.71 x 10⁴ ) を通り、( -0.988, 0.152 ) 方向を向いた直線であり、

と表されます。

y ( または x ) だけの誤差を考えるか、それとも x, y 両方の誤差を考えるかによって、二通りのモデルが考えられました。前者は、x を入力データとみなして、その出力結果 y に誤差があるとして、その誤差を最小にするときの直線を、後者は x, y 両方とも測定データであり、それぞれに誤差があるとしてそれを最小にするときの直線をそれぞれ求めていることになります。どちらを採用するかは、データをどのように解釈するかに依存することになり、どちらが正しいモデルなのかということは言えません。前例において、人口に対して交通事故発生件数が一意に決まるという前提で考えて、測定の誤差(例えば記録されない事故の存在など)によってバラツキが発生すると考えれば前者の方法が、人口にも交通事故発生件数にも同等のバラツキがあって、そのバラツキを最小限にするような近似直線を求めたい場合には後者の方法が利用されることになります。

相関係数には別の見方をすることもできます。N 個からなる x の値によるベクトルを x, 同じく y の値によるベクトルを y としたとき、共分散 s_xy はこの二つのベクトルの内積 ( x - μ_x, y - μ_y ) を N で割ったもので表されます。また、分散 v_x = ( x - μ_x, x - μ_x ) / N = || x - μ_x ||² / N, v_x = || y - μ_y ||² / N なので、標本相関係数 r は

となります。「コーシー = シュワルツの不等式(Cauchy-Schwarz Inequality)」より、r は -1 から 1 までの値をとります。x - μ_x / || x - μ_x || は x - μ_x をその大きさで割ったことになるので、単位ベクトルとなり、y - μ_y についても同様に考えられるので、この単位ベクトルをそれぞれ e_x, e_y とすれば、r = ( e_x, e_y ) と表されます。この二つの単位ベクトルがなす角度を θ としたとき、

なので、r とは x と y のなす角度を表すことになります。

N 次元ベクトル空間上にあるベクトル x, y に対し、それらが同じ方向か全く逆の方向を向いているとき、|r| は 1 になります。また、直交していれば r = 0 になります。同じ方向か全く逆の方向を向いているということは、任意の定数 a, b (但し a ≠ 0) を使って y = ax + b と表されることになるので、xy 平面上にプロットすれば直線上に並ぶことになります。また、例えば

で表されるとき、( Kcosα, Ksinα ) と ( Kcos(-α), Ksin(-α) ) どうしが互いに打ち消し合って、r = 0 になります。K の値をいろいろと変化させれば、平均を原点に、周囲に散らばったような分布を作ることができるので、そのような場合は r がゼロに近い値になります。

ところで、上記の r = 0 になる例では、x² + y² = K² という関係式が成り立ちます。つまり、x と y の間には、線形ではないものの相関があるということになります。この例から分かるように、x と y が互いに独立な場合、

なので相関係数は必ずゼロになりますが、相関係数がゼロでも x と y が互いに独立であるとは限りません。他の例として、y = Kx に従うデータと y = -Kx に従うデータが混在していたような場合、標本を元に相関係数を計算するとほぼゼロに近い値が得られますが、実際には相関のある二つのデータが混在しているものと見なすことができます。よって、相関係数だけで判断するのではなく、相関図で分布の確認をするということも必要であるということになります。

x の確率密度関数が一様分布 p_-1,1(x) に、また、y の確率密度関数が同じく一様分布 p_-2,2(y) に従うと仮定します。ここで、一様分布 p_a,b(x) は

で表されます。x は -1 ≤ x ≤ 1 の区間、y は -2 ≤ y ≤ 2 の区間で等しい確率で分布することになるわけです。x と y が独立であれば、ある x に対して y は -2 ≤ y ≤ 2 の区間で等しい確率で発生します。よって、x と y をプロットしていくと、矩形 ( -1, -2 ) - ( 1, 2 ) 内に等しい確率で分布することになります。このとき、標本相関係数はゼロに近い値になります。しかし、y が x に対して一意に決まるとき、y は -2 ≤ y ≤ 2 の中のある定数を常に示すことになります。測定誤差などによってある程度揺らぐことがあったとしても、それは定数の付近に密集するので、分布は矩形ではなく、ある帯状の形になるでしょう。その形状が直線に近ければ、標本相関係数は 1 に近くなります。

その他に、相関係数の絶対値は x, y が一次変換されても変化しないという特徴もあります。すなわち、

と変換したとき、x, y, u, v の母平均をそれぞれ μ_x, μ_y, μ_u, μ_v、分散をそれぞれ v_x, v_y, v_u, v_v、x,y の共分散を γ_xy、u,v の共分散を γ_xy とすれば、

という関係式が得られます。よって、x,y の相関係数を ρ_xy、u,v の相関係数を ρ_uv としたとき、

となります。

二変数間の相関は、常に線形関係で表されるわけではありません。相関図での例で示した、65 歳以上の人口が全人口に対して占める割合の人口との相関は負の相関がありますが、その関係は直線では示されないように見えます。回帰曲線を使う場合、一次式ではなく二次以上の多項式を利用したり、もっと一般的な関数の和で表すようなことも考えられますが、x, y のいずれかを変数変換することによって、線形関係で表すことができる場合があります。このときよく利用されるのが、対数やべき乗への変換です。

上のグラフはどちらも、65 歳以上の人口が全人口に対して占める割合の人口との相関を示したものですが、右側は x 軸(人口)を底を 10 とする対数に変換しています。両者を比較すると、右側の対数で変換した結果のほうが、線形に近似されている様子が分かります。両者の相関係数は、左側が -0.611、右側が -0.668 となって、対数に変換したほうが相関は強くなります。

標本共分散を計算するためのサンプル・プログラムを以下に示します。

/*
  sampleCovariance : 標本共分散を求める

  const vector<T> &x, &y : データ

  データの個数が一致しない場合はゼロを返す
  データの個数がゼロならゼロを返す

  戻り値 : 標本共分散
*/
template<class T> T sampleCovariance( const vector<T>& x, const vector<T>& y )
{
  unsigned int n = x.size();
  if ( n != y.size() ) return( 0 );
  if ( n == 0 ) return( 0 );

  // 標本平均の算出
  T mx = sampleAverage( x );
  T my = sampleAverage( y );

  // ( x - mx )( y - my ) を求める
  vector<T> xy( n );
  for ( unsigned int i = 0 ; i < n ; ++i )
    xy[i] = ( x[i] - mx ) * ( y[i] - my );

  return( sampleAverage( xy ) );
}
template double sampleCovariance( const vector<double>& x, const vector<double>& y );
template long double sampleCovariance( const vector<long double>& x, const vector<long double>& y );

共分散を求めることができれば、x, y の標本標準偏差で割れば相関係数を、また x, y の分散と標本平均を使って、回帰直線の傾きと切片を求めることができます。

3) 順位相関係数(Rank Correlation Coefficient)

今まで、相関係数として定義していた式

は、「ピアソンの積率相関係数(Pearson Product-moment Correlation Coefficient)」と呼ばれ、測定などによる誤差が正規分布に従うことを前提にしています。また、二変数間の関係は線形であるという条件も必要です。これらの前提条件が成り立たない場合でも二つのデータの相関を表すことができるパラメータとして「順位相関係数(Rank Correlation Coefficient)」というものがあります。

順位相関係数の考え方は単純で、量的データの代わりに順序尺度としての順位を用いて相関係数を計算するというものです。二つのデータ ( x, y ) において、x, y それぞれについて大小関係から順位を付け、得られた順位から相関係数を計算します。人口と交通事故発生件数のデータを例に順位相関係数を計算すると、

より、人口と交通事故発生件数の平均 m_r(x), m_r(y)、分散 v_r(x), v_r(y)、共分散 s_r(x,y) を求めると、

となるので、順位相関係数は 171 / 184 = 0.927 になります。実際のデータから求められた分散 v_x, v_y と共分散 s_xy はそれぞれ 6.81 x 10⁶, 2.63 x 10⁸, 4.05 x 10⁷ なので、相関係数は 0.956 となり、ほぼ同じ値が得られることになります。この順位相関係数は「スピアマンの順位相関係数(Spearman's Rank Correlation Coefficient)」と呼ばれ、イギリスの心理学者「チャールズ・エドワード・スピアマン(Charles Edward Spearman)」によって提唱されました。彼は「因子分析」の先駆者としても知られています。この値は通常、ρ を使って表されます。

データ数が N の場合、順位は 1 から N までとなるので、その平均は m_r = Σ_i{1→N}( i ) / N = ( N + 1 ) / 2 で求められます。また、自然数の二乗和は ( 2N³ + 3N² + N ) / 6 (求め方については「数値積分法 -1-」の「ベルヌーイ数」参照) となるので、分散は v_r = ( 2N³ + 3N² + N ) / 6N - { ( N + 1 ) / 2 }² = ( N² - 1 ) / 12 になります。共分散は、各データの順位を ( r_xi, r_yi ) とすれば

となります。途中で Σ_i{1→N}( r_xi ) = Σ_i{1→N}( r_yi ) = N・m_r であることを利用しています。ところで、

より ( r_xi - r_yi )² = d_i² とすれば Σ_i{1→N}( r_xir_yi ) / N = ( 2N² + 3N + 1 ) / 6 - Σ_i{1→N}( d_i² ) / 2N になります。よって、

なので、順位相関係数の式は

になります。Σ_i{1→N}( d_i² ) は ( x, y ) のペアに対する順位差の二乗和を意味するので、この値を利用することで比較的簡単な計算で求めることができることになります。データそのものの代わりに順位を使って計算しているだけで、相関係数の式としてはピアソンの積率相関係数と変わらないので、必ず -1 ≤ ρ ≤ 1 の範囲の値を示すなど同様の特徴を持っています。

順位を決めるとき、等しいデータがある場合は平均順位を計算します。i 番目から j 番目の順位のデータが全て等しい場合 ( j > i )、i から j までの和をそのデータ数で割れば平均順位が求められ、i から j までの和は Σ_k{i→j}( k ) = { j( j + 1 ) - i( i - 1 ) } / 2、データ数は j - i + 1 になるので、平均順位は

を計算することで求めることができます。例えば、3 番目から 5 番目が等しいデータである場合は

になります。i 番目と j 番目のちょうど中央にあたる場所が平均順位となると考えれば分かりやすいと思います。

順位相関係数の中でよく利用される他の種類として「ケンドールの順位相関係数(Kendall Tau Rank Correlation Coefficient)」というものがあります。"Tau"という言葉が含まれていることから想像できるように、一般的には τ で表されます。この場合も、各データについて順位を求めるところまではスピアマンの順位相関係数と同じです。すると、各データのペア ( x_i, y_i ) ( i = 1, 2, ... N ) に対してその順位 ( r_xi, r_yi ) が得られます。次に、この順位を各データのペアの全組み合わせに対して r_xi, r_yi 両方に対して比較します。比較するデータ列の順位を ( r_xi, r_yi ), ( r_xj, r_yj ) としたとき、r_xi > r_xj かつ r_yi > r_yj または r_xi < r_xj かつ r_yi < r_yj であれば P に、r_xi > r_xj かつ r_yi < r_yj または r_xi < r_xj かつ r_yi > r_yj であれば Q に 1 を加え、r_xi = r_xj または r_yi = r_yj である場合は無視します( P, Q のどちらにも加算しません)。簡単に言えば、順位の大小関係が x, y のどちらも一致していれば P を、逆転していれば Q をそれぞれ加算することになります。x を昇順に並べたとき、y が昇順に並んでいれば、P の値が最大、Q の値が最小(ゼロ)となります。逆に、y が降順に並んだ場合、Q の値が最大、P の値が最小になります。よって、P と Q の差の絶対値は、( x, y ) に相関があるような場合に大きくなる傾向にあります。

具体的な例を示します。国語と数学の試験で各生徒が以下のような順位になったとき、二教科のあいだで相関があるかを調べたいとします。

A から E までの生徒それぞれの組み合わせについて、国語と数学の順位の大小関係が同じであれば P、逆であれば Q とします。例えば、A, B 両者を比較したとき、国語は B の方が順位が高く、数学は A の方が順位が高いので逆の関係となっています。よって、このときは Q で示されます。

P, Q の数を数えるとそれぞれ 6, 4 になります。従って、P - Q = 6 - 4 = 2 と求めることができます。

N 個のデータから 2 個のデータを選ぶ組み合わせは _NC₂ = N( N - 1 ) / 2 になります。同順位がなければ P + Q = _NC₂ であり、0 ≤ P ≤ _NC₂, 0 ≤ Q ≤ _NC₂ なので、P - Q / _NC₂ は -1 から 1 の間の値をとります。そこで、

とすれば、τ は -1 と 1 の間の値をとり、相関係数として扱うことができるようになります。これがケンドールの順位相関係数になります。上の例では τ = ( 2 x 2 ) / ( 5 * 4 ) = 0.2 と求められます。

x の順位は配列のインデックスを表しているとすれば、y の順位が逆転している場合、配列の要素の大小関係が逆転している(ソートされていない)ことを意味するので、Q はその個数を表すことになります。2Q / N( N - 1 ) はその比率を表し、この値がゼロなら全ての要素はソートされていることを、1 なら要素は逆順にソートされていることを意味します。この値のことを "Kendall Tau Distance" または "Bubble-sort Distance" といいます。後者の名前は、この値がソートをするときの交換回数を表すために付けられたようです。

順位相関係数を求めるためのサンプル・プログラムを以下に示します。

/*
  getRank : データ列から順位を求める

  const vector<T>& data : データ列
  vector<double>& rank : 順位を登録する配列
*/
template<class T> void getRank( const vector<T>& data, vector<double>& rank )
{
  // データ型の宣言
  typedef typename std::pair<T,unsigned int> pair;
  typedef typename std::multimap<T,unsigned int> mmap;
  typedef typename mmap::const_iterator mmap_cit;
  typedef typename mmap::iterator mmap_it;
  typedef typename mmap::difference_type mmap_diff;

  unsigned int sz = data.size(); // データ列のサイズ
  rank.resize( sz );
  if ( sz == 0 ) return;

  /*
    データのソート
    multimapへの登録によってキー順にソートされる
  */
  mmap orderedData;
  for ( unsigned int i = 0 ; i < sz ; ++i )
    orderedData.insert( pair( data[i], i ) );

  // 順位への変換
  mmap_cit cit = orderedData.begin();
  unsigned int r = 0; // 現在の順位
  while ( cit != orderedData.end() ) {
    mmap_cit ub = orderedData.upper_bound( cit->first ); // キーが等しい範囲の末尾
    mmap_diff diff = std::distance( cit, ub ); // キーが等しい要素の個数
    // 平均順位の計算
    double rankAvg = ( diff == 1 ) ?
      r :
      ( (double)( r * 2 ) + (double)( diff - 1 ) ) / 2.0;
    // 平均順位を対象位置に登録
    while ( cit != ub ) {
      rank[cit->second] = rankAvg;
      ++cit;
    }
    r += diff;
  }
}
template void getRank( const vector<int>& data, vector<double>& rank );
template void getRank( const vector<unsigned int>& data, vector<double>& rank );
template void getRank( const vector<double>& data, vector<double>& rank );

/*
  spearmanRankCorrCoef : スピアマンの順位相関係数を求める

  const vector<T> &x, &y : データ列

  戻り値 : スピアマンの順位相関係数
*/
double spearmanRankCorrCoef( const vector<double>& x, const vector<double>& y )
{
  // x, y の順位を登録する配列
  vector<double> rank_x;
  vector<double> rank_y;

  unsigned int n = x.size(); // データ列のサイズ
  if ( n != y.size() ) {
    cout << "The size of data ( x, y ) must be the same size." << endl;
    return( NAN );
  }
  if ( n < 2 ) {
    cout << "The size of data must be greater than one." << endl;
    return( NAN );
  }

  // 順位の計算
  getRank( x, rank_x );
  getRank( y, rank_y );

  // 順位差の二乗を求める
  vector<double> d( n );
  for ( unsigned int i = 0 ; i < n ; ++i )
    d[i] = pow( rank_x[i] - rank_y[i], 2 );

  double d_sum = sum( d ); // 順位差の二乗和
  double denom = (double)n * ( pow( n, 2 ) - 1.0 ); // 分母( N( N*N - 1 ) )

  return( ( denom - 6.0 * d_sum ) / denom );
}

/*
  kendallRankCorrCoef : ケンドールの順位相関係数を求める

  const vector<T> &x, &y : データ列

  戻り値 : ケンドールの順位相関係数
*/
double kendallRankCorrCoef( const vector<double>& x, const vector<double>& y )
{
  // x, y の順位を登録する配列
  vector<double> rank_x;
  vector<double> rank_y;

  unsigned int n = x.size(); // データ列のサイズ
  if ( n != y.size() ) {
    cout << "The size of data ( x, y ) must be the same size." << endl;
    return( NAN );
  }
  if ( n < 2 ) {
    cout << "The size of data must be greater than one." << endl;
    return( NAN );
  }

  // 順位の計算
  getRank( x, rank_x );
  getRank( y, rank_y );

  int p = 0; // 順位の大小関係が x と y で同じ数
  int q = 0; // 順位の大小関係が x と y で異なる数
  for ( unsigned int i = 0 ; i < n - 1 ; ++i ) {
    for ( unsigned int j = i + 1 ; j < n ; ++j ) {
      if ( rank_x[i] == rank_x[j] || rank_y[i] == rank_y[j] )
        continue;
      if ( ( rank_x[i] > rank_x[j] && rank_y[i] > rank_y[j] ) ||
           ( rank_x[i] < rank_x[j] && rank_y[i] < rank_y[j] ) ) {
        ++p;
      } else {
        ++q;
      }
    }
  }

  return( (double)( p - q ) * 2.0 / (double)( n * ( n - 1 ) ) );
}

getRank はデータ列からその順位を求めて返す関数です。データ列を引数 data として渡せば、同じく引数として渡された配列 rank に各要素に対する順位を格納して返すようになっています。例えば、data = { 0.5, 0.1, 0.3, 0.4, 0.2 } ならば、rank = { 4, 0, 2, 3, 1 } として返されます。データのソートには STL(Standard Template Library)で提供されている連想配列 multimap を利用しています。キーが一意になる map に対し、multimap は重複を許す連想配列です。実データをキー、配列 data のインデックスを値としてメンバ関数 insert を使って挿入すると、その要素はキーの持つ大小関係によって自動的にソートされていきます。よって、全てのデータを登録した段階でソートは完了した状態になり、反復子(iterator)を使って先頭から順に配列のインデックスを取り出して、その位置に順位を書き込むことで求める順位の配列が得られます。
重複した値がある場合、平均順位を計算して rank に格納します。そのため、現在の反復子にあるキーに等しい範囲の末尾をメンバ関数 upper_bound で、キーの等しい要素の数をメンバ関数 distance でそれぞれ求め、重複要素があれば平均順位を計算して各位置に格納しています。

spearmanRankCorrCoef は「スピアマンの順位相関係数」を求めるための関数です。getRank を使ってデータ列 x, y それぞれの順位を求めたら、順位差の二乗和を求めてから計算式に従い値を求めるだけの単純な処理です。kendallRankCorrCoef は「ケンドールの順位相関係数」を求める関数で、こちらも最初にデータ列の順位を計算し、順位の大小関係をチェックしながら変数 p, q にその数を計算した上で、計算式に従い順位相関係数を求めます。

補足1) ラグランジュの未定乗数法(Lagrange Multiplier)

条件 x + y = 1 の元で、関数 f( x, y ) = 2x² + 3y² を最小にするときの ( x, y ) と、その最小値を求めるというような場合、通常は y = 1 - x を f( x, y ) に代入して y を消去して極値を求めるという方法を利用しますが、変数の消去を使わずにこのような問題を解く手法として「ラグランジュの未定乗数法」があります。

f( x, y ) = 2x² + 3y² は、( x, y ) = ( 0, 0 ) でゼロになり、xy 平面に接した状態になります。また、その形状は楕円型で、ちょうどお椀のように見え、xy 平面に平行な面で切ったときの切り口は楕円の形になります。xy 平面上に直線 x + y = 1 を描画して、その線上に、xy 平面に垂直な平面を描画したとき、その面で切られた楕円型の切り口は放物線となるので、その極小値となるところが最小値となります。ちょうど、楕円形の粘土模型があって、直線 x + y = 1 上にあるひもを上に持ち上げると、粘土が切れたときの切り口が放物線となる様子を思い浮かべると分かりやすいと思います。すると、そのひもが粘土に接したときの点が、求める ( x, y ) とその最小値になります。もし、ひもが粘土に接しておらず、交差した(粘土にひもが食い込んだような)状態になっているとすれば、その点から切り口をたどると値が大きく、または小さく変化するので、その点は最小値とならないからです。

楕円型 f( x, y ) が直線 x + y = 1 と接している点を含むように、xy 平面に平行な面で切ったとき、楕円 f( x, y ) = c (定数) に直線 x + y = 1 が接した形状を見ることができます。この時、f( x, y ) = c の接点での法線ベクトルと直線 x + y = 1 の法線ベクトルは平行になるので、f( x, y ) = c の法線ベクトルを ∇f、直線 x + y = 1 の法線ベクトルを ∇g としたとき、ある定数 λ があって

と表すことができます。f( x, y ) = c と x + y = 1 の接点における、f( x, y ) = c の接線の傾きを求めるため、f( x, y ) = c を x で微分すると、

より dy / dx = -2x / 3y となるので、接線に垂直なベクトルは ( 2x, 3y ) となって、

より、

となります。 x + y = 1 に ( x, y ) = ( λ / 2, λ / 3 ) を代入して λ = 6 / 5 となり、これより

という解が得られ、最小値は 2 x ( 3 / 5 )² + 3 x ( 2 / 5 )² = 6 / 5 になります。

一般的に、制約条件 g( x₁, x₂, ... x_N ) = g(x) = 0 もとでの N 変数関数 f( x₁, x₂, ... x_N ) = f(x) の極値を求める問題を考えた場合、g(x) = 0 は N 次元空間上の曲面を表し、f(x) が極値をとる点において、その点を含む曲面 f((x) = c (定数) が g(x) = 0 と接しています。よって、接点における両者の法線ベクトルは平行になり、f((x) = c の法線ベクトルを ∇f、g(x) = 0 の法線ベクトルを ∇g としたとき、ある定数 λ があって

と表すことができます。このときの定数 λ を「ラグランジュ乗数(Lagrange Multiplier)」といいます。

∇f, ∇g は、偏微分係数 ∂f / ∂x_i, ∂g / ∂x_i ( i = 1, 2, ... N ) を使って次のように表すことができます。

よって、各成分 x_i に対して

が成り立つ x_i と λ を求めればよいことになりますが、未知数が N + 1 個なのに対し、式は N 個で一つ足りないので、もうひとつ制約条件

を方程式に加える必要があります。そこで、

として、N + 1 個の方程式を

と定義すると、F( x, λ ) を x_i で偏微分すれば

なので、左辺がゼロということは ∂f / ∂x_i = λ∂g / ∂x_i を意味し、これが i = 1, 2, ... N で成り立つので ∇f = λ∇g であることを意味します。また、λ で偏微分すれば

となり、左辺がゼロならば制約条件 g(x) = 0 を意味するので、N + 1 個の連立方程式を解くことで解が得られることになります。

直線への当てはめの中で

が最小になるような ( μ_xj, μ_yj ) を制約条件

の下で求めるような場合、ラグランジュの未定乗数法を利用することができます。まず、

として、三つの連立方程式 ∂F / ∂μ_xj = 0, ∂F / ∂μ_yj = 0, ∂F / ∂λ_j = 0 を求めると、

上側の二式から μ_xj = λ_ja + x_j, μ_yj = λ_jb + y_j を三式目に代入して

より

a² + b² = 1 と正規化すれば λ_j = -( ax_j + by_j + c ) なので、

となって、変数を消去して求めた場合と同様の結果が得られます。

前に戻る

タイトルに戻る