JP2007328744A - 回帰分析方法及びその装置 - Google Patents

Abstract

【課題】 シフト項を含む単回帰モデル、更には一般の重回帰モデルに対する評価関数及びパラメタ推定方法を求める。
【解決手段】 多変量データと目的変数との組み合わせからなる行列の分散共分散行列を求め（Ｓ１０２）、その分散共分散行列を固有値分解することによって、複数の固有値及び対応する固有ベクトルを計算する（Ｓ１０３）。そして、これら複数の固有値より、絶対値の大きい順に累積寄与率を計算し（Ｓ１０４）、所定値を超える累積寄与率に対応する固有値及び固有ベクトルから回帰係数を算出する（Ｓ１０６）。
【選択図】 図１

Description

本発明は、多変量データから目的変数を求める技術に関する。

＜記号の説明＞
以下の記述で用いられる記号の定義を与える。

１．Ｄ次元サンプルベクトルｘ，ｙ（英小文字）の要素をそれぞれｘ(1)，ｘ(2)，…，ｘ(D)及びｙ(1)，ｙ(2)，…，ｙ(D)等と記述する。

２．Ｎ個のサンプルベクトルｘ，ｙのサンプル平均値を＜ｘ＞，＜ｙ＞と記述する。

j=N j=N j=N
＜ｘ＞：＝（1/N）Σｘj(1)，Σｘj(2)，…，Σｘj(D)）^T …式（４）
j=1 j=1 j=1
３．分散共分散：
サンプルベクトルｘjの標準偏差Ｓxxは、式（５）で表される。

j=N
Ｓxx＝（1/N）Σ||ｘj−＜ｘ＞||² …式（５）
j=1
またサンプルベクトルｙjの標準偏差Ｓyyは、式（６）で表される。

j=N
Ｓyy＝（1/N）Σ||ｙj−＜ｙ＞||² …式（６）
j=1
ここで、Ｘ，ｙの共分散は、式（７）で表される。

j=N
Ｓxy＝（1/N）Σ||ｘj−＜ｘ＞|| ||ｙj−＜ｙ＞|| …式（７）
j=1

未知の系に対する入力ベクトルｘと出力ｙの組、或いは変量ｘ，ｙの測定値がそれぞれＮ個与えられたとき、これらｘとｙとの間に線形関係が存在すると仮定すると、その関係は以下の式で表される。

ｙ＝θ1^Tｘ＋θ2 …式（８）
この式（８）において、パラメタθ1及びθ2を求めることにより、これらｘとｙとの関係が得られる。ここで、これらパラメタθ1及びθ2を推定する技術は回帰分析技術として知られている。

このような回帰分析の応用例として、特許文献１には、プロセスの状態を推定する方法として重回帰分析が記載されている。また特許文献２では、測色機による分光反射率の測定データから微小見本の分光反射率を推定するために重回帰分析が用いられている。また特許文献３では、積雪量測定装置から得られる積雪情報に基づいて、積雪量測定装置が設置されていない区域の積雪量を高精度に推定するために重回帰分析を用いている。また特許文献４では、３原色を種々組み合わせて作成したサンプルから実測される１次色濃度と再現色濃度とを基に、重回帰分析を用いて再現色濃度を推定している。更に、特許文献５では、走行時の車両状態から路面摩擦係数を算出するために、ブレーキ圧、車輪加速度及び車輪スリップ率から重回帰分析を用いて路面摩擦係数を推定している。

その他、このような応用例は様々な分野において多数存在するが、これらに共通することは、入力ベクトル（説明変数）から出力値（目的変数）への関数のパラメタを精度良く推定する手段として用いられている重回帰分析が主要な役割を果たしているということである。

この推定に使用される数式は、例えば非特許文献１の第１６５ページに記載されているように、以下のように表すことができる。まず式（８）を変形する。

ｙ＝θ^TＸ …式（９）
但し、
｜ 1, 1,…, 1 ｜
｜ｘ1(1),ｘ2(1), ...,ｘN(1)｜
｜ｘ1(2),ｘ2(2), ...,ｘN(2)｜
Ｘ＝ ｜ ... ｜ …式（１０）
｜ ... ｜
｜ｘ1(D),ｘ2(D), ...,ｘN(D)｜
θ＝（θ2，θ1(1)，θ1(2)，…，θ1(D)）^T …式（１１）
この式（１１）のパラメタθに対する推定式は次のように与えられる。

θ＝（ＸＸ^T）^-1ＸＹ^T …式（１２）
しかし、サンプルベクトルの２つの成分の相関が強い場合は、行列ＸＸ^Tが特異に近くなり、式（１２）によって得られるパラメタベクトルの精度が劣化することが知られている。

このような劣化を防止する技術として、非特許文献２及び非特許文献３には、主成分回帰分析法（ＰＣＲ）及び部分最小２乗法（ＰＬＳ）が開示されている。これらはそれぞれ、主成分の分析を利用し、またＸの分散を最大化するという評価基準で部分空間の基底を選択し、それを基に回帰分析を行うものである。具体的には、パラメタの推定式は次のように表される。

θ＝ＶＳ^-1Ｕ^TＹ^T …式（１３）
上式（１３）において、Ｖ，Ｓ及びＵは、Ｘの特異値分解により得られる行列である。

Ｘ＝ＵＳＶ^T …式（１４）
いまサンプルベクトルを構成する各成分間の相関が弱ければ上述の式（１２）を用いて回帰パラメタを推定できる。しかし、相関の強い複数の成分が存在する場合は、上記式（１３）を用いて回帰パラメタを推定することができる。

しかし、最小２乗法に基づくパラメタ推定方法で得られた推定値には一致性がない。即ち、サンプル数をいくら多くしても、真の値との誤差が残ることが知られている。従って、式（１２）或いは式（１３）で推定されたパラメタは、サンプル数をいくら増やしても真の値との間で誤差が残るという問題がある。

一方、非特許文献４は、式（１５）でモデル化される単回帰分析に対するパラメタ推定のための評価関数を式（１６）で表すことによって、新たなパラメタ推定方法を開示している。

ｙ＝θ1ｘ …式（１５）
Ｌ＝（Ｙ−θ1^TＸ）（Ｙ−θ1^TＸ）^T／（1＋θ1²） …式（１６）
この式（１６）をパラメタθ1で偏微分することにより、以下の式（１７）で示すような、パラメタ推定式が得られる。

θ1＝｛−（Ｓxx−Ｓyy）±（（Ｓxx−Ｓyy）²＋４Ｓxy）^1/2｝／（２Ｓxy） …式（１７）
パラメタの推定値は、式（１７）の２次方程式の解として得られる２つのうち適切な方、例えば推定誤差が小さい方を選択すれば良い。このパラメタ推定方法で得られた推定値は一致性を持つ。
特開平６−１１０５０４号公報 特開平６−１１７９３２号公報 特開平６−３０１６６９号公報 特開平６−３５０８４３号公報 特開平７−０１７３４６号公報 S. K. Kachigan ：「Multivariate Statistical Analysis」, Radius, ( 1991 ). W. Wu, R. Manne ：「Fast regression methods in a Lanczos ( or PLS-1 ) basis. Theory and applications」, Chemometrics and intelligent laboratory systems, 51, pp.145-161 ( 2000 ). R. Ergon ：「Informative PLS score-loading plots for process understanding and monitoring」, Journal of Process Control, 14, pp.889-897 (2004). 甘利、川鍋：「線形関係の推定−最小2乗法は最良であるのか？」応用数理, vol.6, no.2, pp.96-109 (JUN. 1996).

しかしながら上記従来技術には以下に述べるような問題があった。

非特許文献４に開示された、式（１６）の評価関数に基づくパラメタの推定方法は、式（１７）に示す単回帰モデルでシフト項のないものが使用されたに過ぎず、シフト項を含む単回帰モデル、更には一般の重回帰モデルに対する評価関数及びパラメタ推定方法について考慮されていない。

本願発明の目的は上記従来技術の欠点を解消することにある。

上記目的を達成するために本発明の一態様に係る回帰分析方法は以下のような工程を備える。即ち、
多変量データから目的変数を求めるために重回帰分析を行って回帰係数を求める重回帰分析方法であって、
前記多変量データと目的変数との組み合わせからなる行列の分散共分散行列を求める工程と、
前記分散共分散行列を固有値分解することによって複数の固有値及び対応する固有ベクトルを計算するベクトル計算工程と、
前記複数の固有値より絶対値の大きい順に累積寄与率を計算する計算工程と、
所定値を超える前記累積寄与率に対応する前記固有値及び固有ベクトルから回帰係数を算出する算出工程と、を有することを特徴とする。

上記目的を達成するために本発明の一態様に係る回帰分析装置は以下のような構成を備える。即ち、
多変量データから目的変数を求めるために重回帰分析を行って回帰係数を求める重回帰分析装置であって、
前記多変量データと目的変数との組み合わせからなる行列の分散共分散行列を求める手段と、
前記分散共分散行列を固有値分解することによって複数の固有値及び対応する固有ベクトルを計算するベクトル計算手段と、
前記複数の固有値より絶対値の大きい順に累積寄与率を計算する計算手段と、
所定値を超える前記累積寄与率に対応する前記固有値及び固有ベクトルから回帰係数を算出する算出手段と、を有することを特徴とする。

尚、この課題を解決するための手段は、本願発明の特徴の全てを列挙しているものではなく、特許請求の範囲に記載された他の請求項及びそれら特徴群の組み合わせも発明になり得る。

本発明によれば、変数ベクトルの分散共分散行列が正則な場合はもちろん、特異に近い場合でも、従来技術より精度が良く、更にサンプル数の増加に伴って真の値に近づいていく推定値を与えることが可能になった。

以下、添付図面を参照して本発明の好適な実施の形態を詳しく説明する。尚、以下の実施の形態は特許請求の範囲に係る本発明を限定するものでなく、また本実施の形態で説明されている特徴の組み合わせの全てが本発明の解決手段に必須のものとは限らない。

図１は、本発明の実施の形態に係る情報処理装置（コンピュータ機器）により実行される処理の流れを説明するフローチャートである。

図２は、本実施の形態に係る情報処理装置（コンピュータ機器）の機能構成を示すブロック図である。

まず図２を参照して、この情報処理装置（コンピュータ機器）の構成を説明する。

図中、ＣＰＵ２０１は、一次記憶装置２０４に記憶されているプログラムに従って、この装置全体の動作を制御している。表示装置２０２は、液晶やプラズマ、或はＣＲＴ等の表示装置で、処理対象のデータや演算結果、ユーザへのメッセージ、エラーなどの表示を行う。入力装置２０３は、例えばキーボードやポインティングデバイスデバイス（マウス）等を含み、ユーザにより操作されて、この装置にデータを入力したり各種コマンドの入力などに使用される。一次記憶装置２０４はＲＡＭなどのメモリを含み、アプリケーションプログラムの実行時、ＣＰＵ２０１により実行されるアプリケーションプログラムが二次記憶装置２０５からこの一次記憶装置２０４にロードされて実行される。またこの一次記憶装置２０４にはＯＳもロードされている。更に、一次記憶装置２０４は、ＣＰＵ２０１による制御処理の実行時に、各種データを一時的に保存するためのワークエリアも提供している。２次記憶装置２０５は、予めインストールされたＯＳやアプリケーションプログラムなどを記憶しており、実行時に、それらプログラムが読み出されて一次記憶装置２０４にロードされる。この二次記憶装置２０５は、例えばハードディスク、ＭＯ等を備えている。通信装置２０６は、ＬＡＮやインターネット等と接続し、これらネットワークとの間でデータのやり取りを行う。バスライン２０７は、上述した各部同士を接続するためのバスである。

尚、本実施の形態に係る処理を実行するプログラムは、予め２次記憶装置２０５に格納されており、入力装置２０３或いは通信装置２０６等からのコマンド入力によって一次記憶装置２０４にロードされ、ＣＰＵ２０１の制御の下で実行される。

本実施の形態では、予め設定した線形係数パラメタベクトルθ1とシフトパラメタθ2に基づいて得られるサンプルベクトルｘj及びサンプルｙjの集合｛（ｙj，ｘj）｝jに対して、従来技術及び本実施の形態に係る回帰分析を実施し、これにより得られた推定値を比較検討する。

サンプルベクトルｘjの集合は以下のようにして算出される。

１．一様乱数randに適当な正実数Ｃ0をかけてξjを算出する。

ξj＝Ｃ0 rand …式（１８）
２．上記ξjに、標準正規分布に従う観測ノイズを付加してサンプルベクトルｘjを得る。

ｘj＝ξj＋Ｎ（0，ε²） …式（１９）
３．ξjの線形関数としてζjを算出する。

ζj＝θ1^Tξj＋θ2 …式（２０）
４．サンプルｙjは、ζjに標準正規分布に従う観測ノイズを付加することによって得られる。

ｙj＝ζj＋Ｎ（0，ε²） …式（２１）
５．上記１から４の処理を繰り返すことによって、サンプルベクトルｘjとサンプルｙjの組の集合｛（ｙj，ｘj）｝jを得る。

以下、本実施の形態に係る処理の流れを図１のフローチャートを参照して説明する。

まずステップＳ１０１で、サンプルベクトルｘi＝(ｘi(1)，ｘi(2)，…，ｘi(d))^Tとサンプルｙiとから新たな変数ベクトルｚiを次式（２２）により作成する。

ｚi＝(ｙi，ｘi(1)，ｘi(2)，…，ｘi(d))^T …式（２２）
次にステップＳ１０２で、このサンプルベクトルｚiから、ｚi（i=1，…，d+1）の平均ベクトル＜ｚ＞を減じて得られるベクトル（ｚi−＜ｚ＞）を横に並べて得られる行列Ｚを求める。これは、
Ｚ＝（ｚ1−＜ｚ＞，ｚ2−＜ｚ＞，…，ｚd+1−＜ｚ＞） …式（２３）
で得られる。

これにより、以下の式（２４）に従って分散共分散行列Ｒを作成する。

Ｒ＝ＺＺ^T …式（２４）
次にステップＳ１０３で、分散共分散行列Ｒの固有値λi及び対応する固有ベクトルｑi（i=1，…，d+1）を算出する。

ここで、この固有値λi及び対応する固有ベクトルｑiは、固有値λiの絶対値の大きい順に並べられているとする。

次にステップＳ１０４で、各固有値λiに対応する累積寄与率ηk を、
j=k j=d+1
ηk ＝（Σλj）／（Σλj） …式（４）
j=1 j=1
により算出する。

次にステップＳ１０５では、予め設定された有効累積寄与率ηCを超える累積寄与率をサーチし（これをηjとする）、そのときの固有値λj+1に対応する固有ベクトルｑj+1をζ1(j)に格納する。

次にステップＳ１０６では、
θ1＝（ζ1(2)/ζ1(1)，ζ1(3)/ζ1(1)，…，ζ1(d+1)/ζ1(1))^T …式（１）
よりパラメタθ1を推定する。また、
θ2＝＜ｙ＞−θ1^T＜ｘ＞ …式（２）
よりパラメタθ2を推定する。但し、ここで＜ｙ＞はｙのサンプル平均値である。

図６は、上述したサンプルベクトルｘi、サンプルｙi、変数ベクトルｚi、その平均＜ｚ＞、分散共分散行列Ｒと、その固有値λi、固有ベクトルｑi、累積寄与率ηk等が一次記憶装置２０４に記憶されている状態を説明する図である。

図３（Ａ）〜（Ｈ）及び図４（Ａ）〜（Ｈ）は、本実施の形態に基づく推定結果を説明する図である。図において、横軸がサンプル数、縦軸が各方法による推定値と真の値との誤差のノルムを示している。

図中、３０１及び４０１で示す線は最小２乗法による結果を、一点鎖線３０２及び４０２はＰＣＲ法による結果を示している。また実線３０３及び４０３は、ＰＬＳ法による結果を示している。更に、点線３０４及び４０４は、本実施の形態に係る結果を示している。尚、図３及び図４では、参照番号３０１〜３０４，４０１〜４０４は簡略化のために他の図３及び図４の（Ａ）〜（Ｃ），（Ｅ）〜（Ｈ）では省略しているが、これらの図においても同様な意味を持っている。

図３（Ａ）〜（Ｈ）はサンプルベクトルの次元が５次元の場合を示し、図４（Ａ）〜（Ｈ）は１０次元の場合を示している。各図３及び図４において、（Ａ），（Ｂ），（Ｃ），（Ｄ）は変数ベクトルの分散共分散行列Ｒが正則な場合を示している。また（Ｅ），（Ｆ），（Ｇ），（Ｈ）は、分散共分散行列Ｒが特異に近い場合の結果を示している。また（Ａ）は、シフトパラメタθ2、（Ｂ）は、線形係数パラメタベクトルθ1に対する結果を示している。

図３（Ａ），（Ｃ）及び図４（Ａ），（Ｃ）より、分散共分散行列Ｒが正則な場合のシフトパラメタθ2に対する推定値は、全ての方法で同等の精度が得られていることがわかる。

また図３（Ｂ）及び（Ｄ）及び図４（Ｂ）及び（Ｄ）より、分散共分散行列Ｒが正則な場合の線形係数パラメタベクトルθ1に対する推定値は、点線（３０４，４０４）で示す本実施の形態が最も良いことが分かる。次いで実線（３０３，４０３）で示すＰＬＳ法、最小２乗法（３０１，４０１）、ＰＣＲ法（３０２，４０２）の順となっている。

一方、図３（Ｅ），（Ｇ）と図４（Ｅ），（Ｇ）より、分散共分散行列Ｒが特異に近い場合のシフトパラメタθ2に対する推定値は、全ての方法で同等の精度が得られている。

また図３（Ｆ），（Ｈ）と図４（Ｆ），（Ｈ）より、分散共分散行列Ｒが特異に近い場合の線形係数パラメタベクトルθ1に対する推定値は、本実施の形態の場合が最も良いことが分かる。そして、次いでＰＬＳ法、最小２乗法、ＰＣＲ法の順となっている。

更に、線形係数パラメタベクトルの推定誤差において、本実施の形態に係る推定値（３０４，４０４）のみが、サンプル数の増加に伴って（グラフの右方向）推定誤差を減少させる傾向にあることが確認できる。

ここでは、「LJ. Peterson, RW. Johnson, CJ. Kerk: "Exploring Relationships in Body Dimensions," Journal of Statistics Education, Vol.11, No.2 (2003)で取り扱われている被験者の身長及び体重の回帰分析を行った。

データファイルには男性２４７人と女性２６０人の身体各部の寸法及び体重が記載されている。

図５は、これら男性と女性の身長及び体重の分布例を示す図である。図において，横軸は身長を示し、縦軸は体重を示している。そして○は男性、△は女性のデータを示している。

これらの内、男性２４７人の身長及び体重のデータを用いて線形回帰分析を行い、身長を体重の一次関数として表現するためのパラメタを推定した。

身長＝（パラメタ１（θ1））×体重＋（パラメタ２（θ2）） …式（２５）
そして、この式（２５）を用いて、女性の体重から女性の身長を推定し、正しい値との誤差の絶対値の平均を算出した。表１は、従来方法と本実施の形態に係る方法に対する推定パラメタ値及び誤差の絶対値の平均値を示している。但し、身長と体重は式（２５）で関連付けられているとする。

この表１より、従来技術としての３方式は推定パラメタ１，パラメタ２が全て同じ値である。従って誤差の絶対値の平均も同じ値になっている。

一方、本実施の形態に係る方式を用いると、パラメタ１及びパラメタ２の推定値が従来方式と異なっている。その結果、誤差の絶対値の平均が従来の３つの方式よりも小さくなっていることが確認できる。

このように本実施の形態に係る方式によれば、線形回帰分析における推定の精度が向上できる。

以上説明したように本実施の形態によれば、変数ベクトルの分散共分散行列Ｒが正則な場合はもちろん、特異に近い場合でも、従来技術よりも精度が良くなっている。

更に、サンプル数の増加に伴って真の値に近づいていく推定値を与えることが可能になった。

（他の実施形態）
以上、本発明の実施形態について詳述したが、本発明は、複数の機器から構成されるシステムに適用しても良いし、また一つの機器からなる装置に適用しても良い。

なお、本発明は、前述した実施形態の機能を実現するソフトウェアのプログラムを、システム或いは装置に直接或いは遠隔から供給し、そのシステム或いは装置のコンピュータが該供給されたプログラムを読み出して実行することによっても達成され得る。その場合、プログラムの機能を有していれば、形態は、プログラムである必要はない。

従って、本発明の機能処理をコンピュータで実現するために、該コンピュータにインストールされるプログラムコード自体も本発明を実現するものである。つまり、本発明のクレームでは、本発明の機能処理を実現するためのコンピュータプログラム自体も含まれる。その場合、プログラムの機能を有していれば、オブジェクトコード、インタプリタにより実行されるプログラム、ＯＳに供給するスクリプトデータ等、プログラムの形態を問わない。

プログラムを供給するための記録媒体としては、様々なものが使用できる。例えば、フロッピー（登録商標）ディスク、ハードディスク、光ディスク、光磁気ディスク、ＭＯ、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、磁気テープ、不揮発性のメモリカード、ＲＯＭ、ＤＶＤ（ＤＶＤ−ＲＯＭ，ＤＶＤ−Ｒ）などである。

その他、プログラムの供給方法としては、クライアントコンピュータのブラウザを用いてインターネットのホームページに接続し、該ホームページからハードディスク等の記録媒体にダウンロードすることによっても供給できる。その場合、ダウンロードされるのは、本発明のコンピュータプログラムそのもの、もしくは圧縮され自動インストール機能を含むファイルであってもよい。また、本発明のプログラムを構成するプログラムコードを複数のファイルに分割し、それぞれのファイルを異なるホームページからダウンロードすることによっても実現可能である。つまり、本発明の機能処理をコンピュータで実現するためのプログラムファイルを複数のユーザに対してダウンロードさせるＷＷＷサーバも、本発明のクレームに含まれるものである。

また、本発明のプログラムを暗号化してＣＤ−ＲＯＭ等の記憶媒体に格納してユーザに配布する形態としても良い。その場合、所定の条件をクリアしたユーザに対し、インターネットを介してホームページから暗号化を解く鍵情報をダウンロードさせ、その鍵情報を使用することにより暗号化されたプログラムが実行可能な形式でコンピュータにインストールされるようにする。

また、コンピュータが、読み出したプログラムを実行することによって、前述した実施形態の機能が実現される形態以外の形態でも実現可能である。例えば、そのプログラムの指示に基づき、コンピュータ上で稼動しているＯＳなどが、実際の処理の一部または全部を行ない、その処理によっても前述した実施形態の機能が実現され得る。

更に、記録媒体から読み出されたプログラムが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれるようにしてもよい。この場合、その後で、そのプログラムの指示に基づき、その機能拡張ボードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行ない、その処理によって前述した実施形態の機能が実現される。

本発明の実施の形態に係る情報処理装置（コンピュータ機器）により実行される処理の流れを説明するフローチャートである。 本実施の形態に係る情報処理装置の機能構成を示すブロック図である。 、 本実施の形態に基づく推定結果を説明する図である。 本実施の形態に係る男性と女性の身長及び体重の分布例を示す図である。 実施の形態の計算に用いるサンプル、変数、計算結果などが一次記憶装置に記憶されている状態を説明する図である。

Claims (12)

1. 多変量データから目的変数を求めるために重回帰分析を行って回帰係数を求める重回帰分析方法であって、
   前記多変量データと目的変数との組み合わせからなる行列の分散共分散行列を求める工程と、
   前記分散共分散行列を固有値分解することによって複数の固有値及び対応する固有ベクトルを計算するベクトル計算工程と、
   前記複数の固有値より絶対値の大きい順に累積寄与率を計算する計算工程と、
   所定値を超える前記累積寄与率に対応する前記固有値及び固有ベクトルから回帰係数を算出する算出工程と、
   を有することを特徴とする回帰分析方法。
2. 前記所定値を超える前記累積寄与率に対応する前記固有値及び固有ベクトルの番号を、固有値の絶対値の大きい順にｊ番目とし（ｊ＋１）番目の固有ベクトルをζ1(j+1)としたとき、前記算出工程では、以下の式（１）（２）を用いて前記回帰係数（θ1及びθ2）を算出することを特徴とする請求項１に記載の回帰分析方法。
   ここでθ1＝（ζ1(2)/ζ1(1)，ζ1(3)/ζ1(1)，…，ζ1(d+1)/ζ1(1)）^T
   …式（１）
   θ2＝＜ｙ＞−θ1^T＜ｘ＞ …式（２）
   ここで、＜ｘ＞，＜ｙ＞のそれぞれは、ｘ，ｙのサンプル平均値である。
3. 前記回帰計数は、ｄ次元変数ベクトルｘj＝（ｘj(1)，ｘj(2)，…，ｘj(d)）^Tの集合Ｘ＝｛ｘ1，ｘ2，…，ｘN｝と、１次元変数ｙjの集合Ｙ＝｛ｙ1，ｙ2，…，ｙN｝に対して、ｙ＝θ1^Tｘ＋θ2となる関係を有することを特徴とする請求項１又は２に記載の回帰分析方法。
4. 前記ζ1(1)と、予め設定された正の実数とを比較し、前記ζ1(1)の方が小さい場合に、前記回帰係数の算出が不可能であると判定する工程を更に有することを特徴とする請求項２に記載の回帰分析方法。
5. 前記計算工程は、前記固有値｛λ1，λ2，…，λ(d+1｝及び前記固有ベクトル｛ｑ1，ｑ2，…，ｑd+1}から、以下の式（３）を用いて、前記累積寄与率が計算されることを特徴とする請求項２に記載の回帰分析方法。
   j=k j=d+1
   ηk ：＝（Σλj）／（Σλj） …式（３）
   j=1 j=1
6. 前記算出工程は、前記累積寄与率が前記所定値を超えるｋを求め、当該ｋを前記式（１）（２）のζ1(j)に代入して前記回帰係数を算出することを特徴とする請求項５に記載の回帰分析方法。
7. 多変量データから目的変数を求めるために重回帰分析を行って回帰係数を求める重回帰分析装置であって、
   前記多変量データと目的変数との組み合わせからなる行列の分散共分散行列を求める手段と、
   前記分散共分散行列を固有値分解することによって複数の固有値及び対応する固有ベクトルを計算するベクトル計算手段と、
   前記複数の固有値より絶対値の大きい順に累積寄与率を計算する計算手段と、
   所定値を超える前記累積寄与率に対応する前記固有値及び固有ベクトルから回帰係数を算出する算出手段と、
   を有することを特徴とする回帰分析装置。
8. 前記所定値を超える前記累積寄与率に対応する前記固有値及び固有ベクトルの番号をｊとし（ｊ＋１）番目の固有ベクトルをζ1(j+1)としたとき、前記算出手段では、以下の式（１）（２）を用いて前記回帰係数（θ1及びθ2）を算出することを特徴とする請求項７に記載の回帰分析装置。
   ここでθ1＝（ζ1(2)/ζ1(1)，ζ1(3)/ζ1(1)，…，ζ1(d+1)/ζ1(1)）^T
   …式（１）
   θ2＝＜ｙ＞−θ1^T＜ｘ＞ …式（２）
   ここで、＜ｘ＞，＜ｙ＞のそれぞれは、ｘ，ｙのサンプル平均値である。
9. 前記回帰計数は、ｄ次元変数ベクトルｘj＝（ｘj(1)，ｘj(2)，…，ｘj(d)）^Tの集合Ｘ＝｛ｘ1，ｘ2，…，ｘN｝と、１次元変数ｙjの集合Ｙ＝｛ｙ1，ｙ2，…，ｙN｝に対して、ｙ＝θ1^Tｘ＋θ2となる関係を有することを特徴とする請求項７又は８に記載の回帰分析装置。
10. 前記ζ1(1)と、予め設定された正の実数とを比較し、前記ζ1(1)の方が小さい場合に、前記回帰係数の算出が不可能であると判定する手段を更に有することを特徴とする請求項８に記載の回帰分析装置。
11. 前記計算手段は、前記固有値｛λ1，λ2，…，λ(d+1｝及び前記固有ベクトル｛ｑ1，ｑ2，…，ｑd+1}から、以下の式（３）を用いて前記累積寄与率が計算されることを特徴とする請求項８に記載の回帰分析装置。
    j=k j=d+1
    ηk ：＝（Σλj）／（Σλj） …式（３）
    j=1 j=1
12. 前記算出手段は、前記累積寄与率が前記所定値を超えるｋを求め、当該ｋを前記式（１）（２）のζ1(j)に代入して前記回帰係数を算出することを特徴とする請求項１１に記載の回帰分析装置。

Images