連載講座
111
川川|川
11 川
11川
11川
11川
11111
川
川
11川
11川|川|川|川
11 川
11川
11川
11111
川
l川|川l川|川|川|川
川
11川
111111川川
11川|川|川|川|川
11 川
11川
11川
11川
11川
11l川川
11川
11川
11川
11川
11川
11川
1111川川川
11川
11川
11川
11川
11川|川
11 川|川
11 川
11川川
11川
11川
11川|川
11 川
11川|川
11 川
11川|川
11 川
11川
11川
11l川|川川
11川|川|川|川|川
11 川
11川
11川
11111
川
川
11川
11川
11川|川|川
11 川
11川|川
11 1川
11川|川
11 川|川11 川|川
11 川
11川|川
11 川
1111111
川
川|川|川|川|川|川|川11
川
11川|川"川"川
11 川
11川
11川|川|川|川|川
11 川
11川"川川"川"川|川11 l川|川|川|川|川|川"川"川"川"1
川"川"川"川|川"川"川|川"川"川"川"川"川"川"川"川"川川"川"川"川"川"川"川"川"川"川11
川"川"川"川"川"川"川"川"川
11
"1
川川"川
11 川"川"川"川"川"川"川
11
川"川
11 川
11川
11川"川"川"川"川"川"川"川11
川"川"川"川"川"川11 川
11川
11""1l行列演算用言語 LAMAX-S(2)本郷茂，八巻直一，内田智史
1.はじめに;文法に対応させる簡易コンパータが用意されていますの

前回では、 LAMAX-S の概要と比較的簡単な LAMAX

で、必要な方に配布できるよう準備しています。もちろ

S プログラムについて紹介しましたが、今回は、引き続ん、プログラムロジックなどは、 LAMAX-S(PC98 版)をき LAMAX-S のプログラムを通して、よりきめ細か L 、文JH いてチェックしであります。

法を説明して行きます。
今回からは、行列表現の数式をそのまま LAMAX-S で占きドすための仕ブi を示すために、本格的な 3 つのプロ

2 ，ヤコビ法による固有値・固有ベクトル
国有値・固有ベクトルは、 LAMAX-S の組込み関数でも求めることができますが、本節では LAMAX-S の例題

グラム例で紹介します。
教育的な配慮のために、効率を考えす叫に行列l の数式をとして、災対称 1 j!i'U の全間有値とそれに対する同台ーベク

そのまま LAMAX-S で表現しています。しかし、我々は、トルを求めるヤコビ法のプログラムを紹介します。このLAMAX-S の最終目標として、高度な最適化機能を実現

プログラムは、メインプログラムと固有値・固有ベクト

し、極力効率を落とさないように動作することを狙ってルを求めるサブルーチン (.JcbMtd) 、およびいくつかのいます。したがって、今回ここで紹介するプログラムはす

補助サブルーチンからなります。メインプログラムでは、べて将米的な LAMAX-S の最適化機能を期待して、特に(1) 単にデータをセットして、 (2) サブルーチン JcbMtd

効率化のための考慮は払わず、記述性を重視しています。

を呼び出して固有値・固有ベクトルを計算し、 (3) その

まず、最初に、固有値・固有ベクトルをヤコビ法で求値を表示しているだけです。実際の計算は、サブルーチ

めるプログラムを紹介します。後半のプログラム例で紹ン JcbMtd で行われ、その内容は、次のヤコビ法のアル介する多変量解析では、固有値を求めることになります

ゴリズムに従っています。 h )il] A の固有値・固有ベクトが、立立初からライブラリを利用するよりは実際に自分でルを求めるとします。

I，'d 干J11!( を求めるアルゴリズムを確認するためのプログラ|ステップ oミングが必要となるでしょう。ここで紹介するプログラ

める。 EVC の初期値として単仇行列をセットする。

ムリストを見ていただければ、 LAMAX-S のプログラム|ステップ 1と実際のアルゴリズムとの一致性を確認していただける

q ，で絶対値最大のものを探す cと liiJI~，Hこ、プログラミングがiィ、慣れな人でも教科書(教科う。次に多変単解析の'1'から、主凶:時分析と主成分分析
を取りあげ、 LA J\<IAX-S での基本的な統計社算の仕方か

ら多変盤解析の具体的なプログラミングを示します。
なお、 L:\:-'lAX-S の文法は、 1992 年の日本 OR 学会

事例研究奨励賞ソフトウェア部門授賞後も、フルセット
LAMAX-S 開発に向けての改訂作業が行われていますっ

ここに IJ~ すプログラムは、全て改訂された LAMAX-S の

来rr 文 l，t に沿ってコーディングされたものであり、そのままでは、 LA t. IAX一 S(PC98 版)では動作しないことをお附りしておきます【実際に正IJ I'I'させるためには、新しい

ほんごう

ここで I が行、 j が列の番号を示す D|ステップ 2 1I胡;角 θ の計算と座標変換行列 P の作成
1,u主S Cl IJ∞ti(} =一三 .11 +::.よSill(}二一一一二一
ゾ2(
3'
2,13 cos(
}

V

ただし，

α= 円 aij/
3=Ja2+arj

ε== sign(αii ‑ ajj)
これを用いて，次の座標変換のための行列 P を作成する。 p は Piゎ lうj 以外の対角部分は1. 0 で、それ以外のしげる専修大学

長京は 0 の行列ですが、 Pii〒 214 川崎市多摩区三田 2-1-1

やまき

なおかず脚システム計画研究所

うちだ

さとし神奈川大学

15
4 (46)

1 t= t+ l として，行列 A の非対角要素のmaxaij(i 千 j)

冷はたくさんありますが、今回のプログラムでは [1 ])を凡ながら容易にコーディングできることが分かるでしょ

1t ニ o として，固有ベクトルの初期点を決=cos(} , Pi
j ==sinO , Jう 2 二

一日 11 0 , Pjj =c刷。となる行列です。

|ステ y プ :1 |行刈 A の兇新と間有ベクトルの更新© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オベレーションズ・リサーチ

A =PAp T 、

EVC=EVC ・ pT

LAMAX-S でのサブルーチンの使い方は、 FORTRANに完全に準拠しています。すなわち、行列(ベクトルを含|ステップ 4 |収束チェック

む)をサブルーチンに引き渡す場合には、その構造、要(1) 行列 A が数値的に十分対角行列になったら正常終了(固有値は A の対角部分)

点の型、数学的性質(正定値など)が完全に一致していな
ければなりません。

(2) 反復回数が n 2 を越えたら終了(失敗)サブルーチン JcbMtd では、変数 A ， EVL , EVC が、J'

|ステップ 5 |ステヅプ 1 へ戻る。

法 II 行 n 列の行列として宣言されています。寸法Jlは、

c

ヤコビ法:固有値・固有ベクトル(メインルーチン〕

親のルーチンから受け取るので、この宣言は FORTRAN

parameter(eps=1e‑6)の終合配列に似ています。

cpar四eter(n=4)<*行列の次元数

parameter(epssub=1e・ 3)<*計算誤差の許容値。一致比較の許容値<*最大繰り返し回数

parameter(maxcnt=20)

real:matrix[n , n:symmetric] A 併求める行列
EVC
制固有ベクトル
real:vector[n] EVL<*固有値
A
I1 , 2 , 3 , 4付固有値・固有ベクトルを
real:matrixln~n]c~) ~) ~)?

3; 6; 8; 9ccc<*求める行列のセット

st~p

で P は、 n にどんな値が代入されても良いように、動的

subroutine JcbMtd( A, EVL , EVC , n , eps , ierr)
己注意(1)行 1 1) A の固有値 固有ベクトルを求める(2) 行夢1I A は，対称行列である.
(3) 行列 A は，破壊される。

cc

ユnte区er n本行列の次元数

real:matri玄[n， n:sy皿etric] A <事求める行列real:vector[n] EVL<牟行斉1I A の固有値
:real: 回trix [n， n) E
VC.
..<*行列 A の国有ベクトルreal ep5<*計算誤差の許容値
integer ierr <*エラーフラグ 0: 正常、 1: 収束せず
real:
m
a
t
r
i
x[
*J
*:tr回sfor凶 P 制座標変換行列
real alpha , beta , es 柑 θ 計算用変数
logical isdiag<市対角行列かどうかチェック
integer t<場繰り返し数

白=======ステッ 7・ 0 初期設定cがプログラムの実行中変化しません。静的行列の寸法は定数式でなくてはならないので、 P の寸法の宣言に n が

en
a
.

ー

と呼び、プログラムの実行中に寸法が変化しますロこれに対し、寸法が明示された行列を静的行列と呼び、、Nl

4 , 7; 9 , 1
0I
call JcbHtd( AA , EVL , EVC , n , eps , ierr)call mprint(EVL) 。固有値の表示
callmprint(EVC) <*固有ベクトルの表示

ヤコビ法:固有値・固有ベクトル{サブルーチン)cドlÆ標変換行3日UP の宣言で、*は P の寸法を意味します C

LA l\ IAXろでは、このように宣言された行列を動的行列i記述できないのです (n は仮日|数です )0 A , EVL , EVC:

は仮ヲ i 数であるので n を寸法として記述できますっそこ
liJlj として宣三されています。動的行列は、データが代人された時点か allocate 文によって明示的に指定された時点で記憶領域が割り当てられます。 frec 文によって明
示的に領域を解放することもできます。
また、 P [*:transformJ の transform は、これが座伎一変換行列であることを意味します。これについては後述します。 P

1::diagonal[n J は、 P に単位行列をセッ卜することを意味していますっ LAMAX-S では、 l 値一構造」という表現を同値要素行列と呼び、これは要素・がす
べて「値 J のその構造をした行列を意味しますったとえ
ば、 0: :matrix[2 , 3J は 2 行 3 JIJ のゼロ行列を意味し

fV2;1 ・ :di 喝。 nal 同

1
0c
o
n
t
i

ます o 1::matrix [1 0 , 10 :uppeLtriJ は 10 1
j10 タ IJi

で上三角部分の要素だけが 1 であとの要素はゼロのわ JIJ

c=====雷同ステップ 1 非対角要素の中で絶対値最大のものを深す。
t
t +1

ã5m

..xr(A , ‑diagonal)

asm日cC A ，‑dij
aonal)

c======== ステッ 7'2 回転角 θ の計算と座標変俊行列 P の作成

alpha= (Aﾇi， i] ・ A[j , j])!2
beta= sq;t(alpha*ホ 2+HI ， J]**2)e
5= sig( 1. 0ëQ , _A[i , i]-A[j , j])P 1:: 孟iagonal[n]P[i , i] = s 吐rt( (1+e5 冷 alpha!beta~ 本 0.5)
P[i ,jJ = (ë~* A[~， j] )/(2*beta事 P[ェ，ュ] )p[ i.ュ・P [i, i]を意味します【\\:;luaxr 関数は、絶対 11立が放火の行列要素の値のある
iJの位置を返します。出nwxr I対数には色今な制限や条例を付けることができます。拙 maxr(A ，-diagonal) !ま「対角1!!~を除く J という条件ですっ aS J1l axc 関数は同様に列P[J, j] = P[i , i]c======== ステソプ 3 行列 A の更新と固有ヘクト J L.の吏新の位置を返します。 asmaxr と耐 maxc が別々に記述され

ç;====::::=::: ス子ッ 7・ 4 収束チェックていても、 LAMAX-S の最適化機能によって最大値の検dv
EP釘t~pJif( isform(A , diagonal , e
p
5
)
)g
ot
o2
0
if(t ・ gt.n*n)go t 30<*~l星回取 h・ 2 なら終 fc======== ステッ 7・ 5 ステップ 1 へ戻る。

go t
o1
0

と--お話:=::ステッ 7・ 6 正常終了

20 continue:l e l:'rυ
EVL= A七 0> 。固有値のセット
ret

c=======ミステァプ 7 異常終了

30 cont ユ nue
lerr = 1
r
e'
t
urn
en<l

1993 年 3 月号

栄は 1 !支しか行われませんので、効率は低下しません。

1凶 fO l'l1l関数は、与えられた行列の非ゼロの配置が特定の構造になっているかどうかを判断します。 isform( A ,山lιollal 、叩 s )は、行列 A が与えられた計算機イプシロ
ンの'1 1 で対角行列になっているかどうかを判断します内
A<O> という表現は、行列 A の対角部分を取りだし、そ© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

(
4
7
)1
5
5

れをベクトルとして扱うことを意味します。 EVL=A<O>

は、行列 A の対角部分をベクトル EVL に代入することを意味します。

このプログラムの特徴は、座標変換行列 P にあります。
教科書の記述と全く同じに記述するという目的でプログ

ラムを作成したので、 P という座標変換行列を作成した訳ですが、実際には行列の回転を行うだけであり、計算量はそれほど多くありません。 p*帥?という計算をその記述どうりまじめに行うと、采算が 2 問必要となってしまいます。しかし、座様変換行列という指定をすることで、 LAMAX-S プリプロセッサの最適化機能によってと記述します。ここでの vector は列ベクト Jレを.

r
v
e
c
t
o
r

1ま行ベクトルを意味しています。

同様に 1 変量の分散は， 11 f聞の 1 を各要素とする n 次元ベクトル ln=(1 , 1,… , l)t の列ベクトルを用いて偏差の2 采和で表現します。

分散 S2=j(x 一言 ln)t(x 一言 ln)
r
e
a
l
:
v
e
c
t
o
r
[
n
]X
Xrn ,S2real
Xrn 1:: 主vec七。r[n]*X/n
S
2 (
X
‑
X
r
n
*
1
::vector[n])' ホ (X-Xrn*1::vector[n])無駄な計算がなくなります。

1::vec 七 or[n] 咋m の記述で全ての要素が平均値とな
3. 統計の初歩の計算と重回帰分析・主成分分析
多変量解析の教科書の解説で代表的な重回帰分析と主る 11 次元列ベクトルが計算でき，その値を 11 次元ベクト
ル x から引くことで偏差になりますの

成分分析を題材にして，それぞれの解析のアルゴリズム

次に，データが n 個で p 変数の多変量の場合の平均とを行列表現のまま LAMAX-S でプログラムしてみましょ分散共分散行列について表してみましょう。多変量の平う。ただし，おことわりしておきますが，この節ではプ

均と分散共分散行列は変量の場合のベクトルの表現

ログラムの記述性を重視しましたので，実行効率や数値が行列に変わりますが. I< .J じ様な記述により去すことが計算上での誤差が少なくなるような工夫はしてありませできますっいま，多変量のデータ X をEEEEEEE'E

、、

ppP

，/

分散 S=i(X-1nJ仰ー l nxp xd)
ただし ， lnx Jl は全ての要素が 1 の(71， xp) 行列で .

、、 EEEEE
F'

一%­

『『』』

AMu‑

U九
M
一

一一

，，a
，a
，
E
・
E・-aaa-‘
-、
-、
-、
z

上式の X-l nxp xd は，各変数の偏差を計算していまは p 次元行ベクトルです。

.!.xtx-xtx
Sニ n

各要素が全て 1 のベクトルを. LAMAX-S では

1
:
:r
v
e
c
t
o
r
[
n
]

平均値ベクトルと分散共分散行列を LAMAX-S で記述してみましょう e

15
6 (48)

xdは対角要素が平均値の (n xp) 対角行列です。

す。次式でも分散共分散行列を求めることもできます。 x

X
r
n r
o
w
s
u
r
n
(
X
)
/
n

1
::
v
e
c
t
o
r[
n
],

12:‑q
z
zZ

.!.ln X副

いて記述することもできます。

勺e---n

平均宮

r
e
a
l
:
v
e
c
t
o
r
[
n
]X
r
e
a
l
X
r
n
X
r
n 1
:
:
r
v
e
c
t
o
r
[
n
]
*
X
/
n
あるいは，行和をもとめる LAMAX-S の組込み関数を用

咽・
4の (n xp) 行列とすると，平均は p 次元ベクトルとなり分散共分散は (p xp) 行列 S となります。

Iz
ｭo
‑

平均言 = ln-xlη

2294

式を LAMAX-S では次のように記述します。

zzZ

行ベクトルとを用いて表現することができます。その数

nの 1 を各要素とする n 次元ベクトル ln =(1 ， 1 ，…， 1) の・ベクトル x ==(X I, X21X3 , …] Xn)t の列ベクトルと. n1聞

12:

散共分散行列や相関係数行列をいかにして LAMAX-S で111，

まず最初に，統計の初歩の計算である 1 変量の平均・分散などの基本的な表現から，多変量を扱う場合での分1 変量での平均は. n 個のデータを成分とする n 次元、
ftsEEEEEz-、-

一一x

きの LAMAX-S での工夫の仕方は別の機会にします。

表現したらよいかについて説明します。

zzzl

ん。実際に，アプリケーションプログラムを作成すると© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オベレーションズ・リサーチ

X

real:matrix[n , p]

で求めます。なお、ここで、 Sij は変数 i ， j の共分散、 Sii ， Sjj

real:matrix[p , p:diagonal] Xmd

は変数 i ， j の分散です。この 1'ij を要素とする (pxp) 相関

real:matrix[p , p]S

係数行列 R は，ゾ訂 1 ，・・・，♂pp を対角要素とする (pxp)(l::rvector[n]*X)'

対角行列 D を用いて，次のように表現するとができます。

Xmd<O>S(
X ‑ l::matrix[n , p]*Xmd)'*

R =D‑1SD‑1,(X-l::ma:七rix [n pJ *Xmd)!nこの式を LAMAX-S で記述してみましょう。

あるいは

DRD

目、，

可晶*
免日明

司J
・V

r

刊、M

+
h
v

，，‘、

、，

時四

qx

・・

冒

3

。

Lgn二

gorkt4CX)

am
可ム白"日司、，，
dsS
1Jwf

rt』

、キ
司

キ

r目、

円、“

ここで. DiagSqrt は分散共分敢行列の対角要素だけを取りだしその平方根を求め結果を列ベクトルとして戻す関数です。この結果は，各変数の標準偏差となっています。

DiagSqrt 関数は. LAMAX-S の組み込み関数ではありません。各自で作成してみてください。この他に.データを規準化(平均 O. 分散1)してから相関係数行列を求めることもできます。具体的な記述については，ここでは省略しますが，実際のプログラム例は主成分分析の最多変量解析の基礎的な数式を行列表現で記述した本も・多く出版されています [2 ， 3] 。また.

FORTRAN,BASIC

一句などのプログラミンクP言語で書かれたプログラムリスト

変数の平均偏差は. X ‑ l::matrix[n ， p]*Xmd で求めています。 1 ::matrix[n , p] は，全ての要素が l の (n X p)

1::matrix[n , p]*Xmd

一引叩一円一z

一Z-h一向:・一UM

Ill

-z 一
z.

一一

\st
JrgEEEEEE、'x¥BEll‑‑/の演算結果は，

1Jrr*

初j の部分に表現してありますので参考にしてください。

ここで.転置の.を付けないと正しく実行されません。各行列を表しています。したがって.

』

0

、l
、 BEl--'sI/

ワM

ハ
Hu

‑Z

一れ

/'glazi‑‑¥

Xmd は，次のような (p xp) 対角行列になります。

nyny*nua(1::rvector[n]*X) ，と記述すると，その演算結果のP ゐ potDと対角行列の対角要素に代入することはできません。 Xmd、

クトルの要素を代入するとき. 1::rvector[n] はで求めた平均値(行ベクトル)を列ベクトルにしてからでない-oa+-T。、 J

LAMAX-S て、注意すべきことは，対角の要素に平均値ベ・-­を説明します。 Xmd は対角行列として宣言しておくと対角以外の要素は全て O とみなして演算してくれます。

xzrg

牟

列に持つ変数 Xmd について LAMAX-S の記述の仕方

キ*

<

VAnu

今までの節でも LAMAX-S の文法を述べていますが，多少解りにくい箇所がありますので，この平均値を対角行

rrtJr‑ Xm'Xmaae

X' キ X/n

τよ可晶〉

rOliJsum(X)/n

-ム

Xm
S

rtLrt

s

s

real:matrix[p , p]ilona1x

real:matrix[p , p]eee,,
ttcXDrreal:matrix[n , p]

Xm

mmvs=

X

aaa=o=
rrrSDR

real:matrix[n , p]real:rvector[p]も一緒に掲載されている本もありますので参考にされるとよいでしょう [3]0

さて，今まで書いてきました平均，分散，相関係数行

声 IJ の LAMAX-S プログラムを使って，重回帰分析と主成分分析プログラムを作ってみます。どちらも.

100 行

程度のプログラムで表現できています。計算方法の詳細についてはここでは述べません。具体的には，参考文献[2 ， 3] などをご覧ください。 LAMAX-S で記述したプロ
グラムは，参考文献問の内容を基に作成してみました。
重回帰分析と主成分分析プログラムの説明は，紙面の結合で詳しくできませんが，プログラム内に注釈を入れの (n X p) 行列となります。

統計の初歩の計算の最後として. p 変数の杵1 関係数行
列を LAMAX-S で表現してみましょう。相関係数は，ておきましたので参考にしてください。

例示したプログラムでは，次のような内容の計算が行われています。

G

T

一一一一~

'
J ‑ ..;s;;..;可71993 年 3 月号

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

(
4
9
)1
5
7

|重回帰分析のプログラム概要|

C 一一

重回帰分析

real*8: 回.trix[* ，ホlreal・ 8: r
vector[*]Xr-データの平均，分散共分散行列.相関係数行列

C日 RR

real*8 ・ matri玄[事.*]国real稼 8:vector[*]real 事 8

R
制重相関係数日<*自由度調整涜み決定係数
町付自由度調整済み重相関慌数

real*8real*8:matr エ玄[*.*]
real・ 8 ・ vector[*J

xty の係数 b を連立」次方程式を

real*3real*8real*8real*8解く solve 文を用いて求めます。

-重相関係数，推定値 b の標準偏差.

t11宣

。残差

56M2 -0 残差の分散
SGK <*残差の標準偏差
VARB <* b の分散
SDB <* b の標準噌差Tv<*β ュ =0 検定の t 値

real 牟 8

XtXb<*基準変数の予測値

YE

real 傘 8:vector[ ・ 1・偏回帰係数の推定値の計算 :lE 規方程式<*基準変数

榊説明変数<*偏回帰係数 {β) の推定値 (b)X

real*8:vector(*]real*8:vector[キ 1&'数 (X) をよさび重凶帰分析<*データの相関係数行列“作業用変数Dyreal*8 ・ vector[*]
real事 8:matrix[* .
*J.データの中から規準変数(目的変数 y) と説明

牟データ

H
E
A
l
l <*データの平均値
SD
榊データの標準備差
COVAR <*データの分散共分散行列

realキ 8:vector[申]real唱団.atrix(* ， .J
real時:田町 ix[* ， *J

Fv

本

目

integer
real*8:vector[・ lintegerβ2=

胃 βF=O 検定の F 値

付決定係数

IS(100) 制作業用変数

創刊!l qrt <*対角要素の平方根を求める関数

N， P~POr~~d(*.*)_)l .p<* N: 観測データの数、 P: 変数の数

allocate Xr:matrix[N.
P]
read(*. 吋 ((Xr [i , j J, ì=l ,P) , i=1 , N)

|主成分分析のプログラム概要|

HEAN

玄O宵 sum(Xr)/dble( 的制平均の計草COVAR = XrJ.Xr/dbleOO ‑ MEAN' 傘 HEAB 榊分散共分散行列の計算
SD= diagsqrtCCOVAR.P)<*標準偏差の計算
O.OdO::matrix[
P.
P]
D<O> = SD
CORR(0** ト1))刊日 VAR* (0**(-1))<*相関係数行列の計算

・データの規準化，相関係数行列

write(*. 吋'データ'
call mp士 int(Xr)-固有値，固有ベクトルwriteC*. 吋'平均値，標準偏差'

call m:p rint け阻 AN 勺 S副〉• EIGENVV プログラムリストは省略しています0write(*. 吋'分散共分散'
callmp i: iD.t(C_Ö~a?

首町玄 ite(収牟円，*吋).相関係数，

首ritel 傘，事)，----ー司自向---------------・・・・胴・・・・・・胴----，

f
品叫1口1 f
ア玄一一t(ωc∞日凹RR
000 COl忍，t エ ewri七 e(* .吋'重回帰分析

write(*.吋' 基準変数の番号
read (*~*)
NDEP
ユf( NOEP .
1e.0) 箆 oto 2000

・主成分負荷量，主成分得点

町 i~eÇ*. り'
町 ite(* 〆)LAMAX-S の文法を含めてより具体的な応月!として行列20プログラムの説明では多少不足しているところもありまていただけたかとおもいます。私共は，この LAMAX一円

'

X
:
;
:
1
.
0
d
O
::1I泡Il~旦手七玄i玄 [H ， P]
〈本 1 カラム目は定数項とする
do 20 i:;:1 , PO
X[1:N.i+1J :; Xr[1:N ，工 S(i)]ホ説明変数の抜き取り
continue
solve (x'*x) キ B = X' 可制偏回帰{高数 β の推定値 b 円計算
YE ~ X*B<*基準変数の予測値の計算
UE :
;Y:
‑ YE
<*残差の計草R2 = (YE' 叫ー〔目時四 (Y l**:j )!dbl.(N) )
制決定係数! ( Y' 叫ー(目間四 (Y)**2)/db1e(H))
R
SQRT <R2)
引童相関係数
回 2= 1. 0dO ‑ db1e(N-1)!db1e(N-P)* <1. O
dO‑R2) <*自由度調整
RF = OSQRT(RF2)。旦由度調整潰み重相関係数
SGM2= 田川 UE ! db1e(N-円。残差の分蔽
SGM = sqrt(SGM2)。残差の槙準偏差VA阻= SGH2!(X' 叫}を用いて教育の現場で優れた教育効果があげられるよう付 b の分散

SDB = diagsqrt(VARB , P)
科 b の標準偏差Tv =B Y. SDB
・ t値
Fv =(R2!dble.ω:1) 2!((!,OdO-R2) !dble(N-P)) 制 F 値使っていただきたいと願っております。次回は，本連載

冒 rite(* ，吋'君事明変数;観測値の数(1)，変数の数【P):' ， N ， P

call mp:dnt(X)の最終回となります。最終回では，数理言 j.ì酉の事例を紹

宵 rite(*.*)，偏回帰係数 β の推定値 b ， b の標準偏差，t 値'

write(*.*) ~.値の自由度【 H~P)=' , N-P
call mp:dnt(
lB.SDB.Tvl)‑ ‑̲.̲‑.
‑.
町 ite(*.*) ，決定係数, , R2
write(*.*) ，重視関係数
・， R

介して，さらに，今後の LAMAX-S の将来構想についても述べたいとおもいます。

P

町田圭

臨H'

冒a

RFZP

zV4

値

差準

FH

t

APHシa
，V
a動

数測
数係め予

係関トの

hho

定相し数
決重」変

度度ののの変涜済度偏差
整整由散準わaa'rnue

調調自分標数阻

FF,

由由値値差基準VH
残残基げ

.,

aa

自白

$e.，， e，*阻0

2

0
0
0
広

{42

他.多変量解析の基礎，東洋経済新報札::;;buwwwwwwwczoC

戸川隼人:詳解数値計算演習、共立出版、 1980

、 f、 r 、r 、r
イif-F

、店
eeeeeeeｭ-ιtttttlot
‑
E
z‑
zz
rl
ri
z‑
z‑
zl
at
oE

参考文献

社.

. ‑ -,

ι..

分析用に取 2.:込んた智明実数の番号ソ， (IS(I) , I=1 , PO)

μ=v“+

すが. LAMAX-S で記述したプログラムの記述性を凡

3. 柳井春夫

ーーー.二

基準変数の番号 ;J .
Il
DEP

N主 [!:H ， HOEP] <*説明変数の抜き取りの固有値解法と多変量解析のプログラムを例示しました門

1
9
7
2

,

writ l!(*.吋'ムを対比しながら例示しましたが，今回の第 2 回目では，

2. 竹内啓PO ，(工 S(工:工=ï;po)-write(*~ 吋'・・・困

第 1 回では. FORATAN と LAMAX-S のプログラ

1

分析用に取り込む劃明変数の個数と番号..

read <*~り4. おわりに

他:多家量解析ノ、ンドフゃック，現代数学

1
9
8
6

1
5
8 (50)© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オベレーションズ・リサーチ

c---- 主成分分析

integern.p

目 .1・8:mat :r i玄白井] Xr

real*8:matrix[*~・] Xp
al*8:matrix[*.*] Xd
e.1・ 8:rvect~r[*]- MEAll:real時 :matr江[*.*] S<*データ数，変数の数

榊データ

real*8:matri玄 r.:・] Sdrea
1*8:matrix[*:*] Cov

借偏差
制基準化
<*変散の空均値<*変数の標準偏差
制作業用変数
制作業用変数<*分散共骨散行列

re a1時 :vector[*j ‑ C
RATE
α皿CR
re a1崎市atrix[*~ 吋 FScore
re a1時:回.trix[仁*] F
Load

<*寄与串
<*累積寄与率<*主成分得点
制主成分負荷量

1・8:matr江[九時 Hd

e a1柑 :m叫口玄ら :*j R

UNIX ワークステーションによる科学技術計算<*相関係数行列

real*8:vector[*]

EVAL
real時 :matrix[*~*] EVEC

出品時 :vector[*]

ハンドヂァク<*坦有値併固有ベクトル制固有値の累積榊固有値の平方根(対角行要1)).1時:皿ã.trix[*: ヰ] CHIDF 制力イ 2 乗値と自由度
玄e.1時 ~;~~t~r[*j diagsqrt <*対角要素の平方根を求める関数
real 時 :vector[*j開ï.8;~~t~i~[* :寸 D
.1 時

間 .rml.rm2 ， r m3.由。作業用変数

本書は，近年のコンビュータ技術の進歩により生み出された低価格・高性能の ws を，十分に活用するため・普通の参考書の 2 倍以上の買数を使って，・最新技術をすぐに役に立つ形で鮮践し，

eC 冨簡によるプログラム例を 80本収録，・そのフロッピ・ディスクを標準添付した，
理工系研究者，技術者，院生に必携の書.
レ

‑
h
r

今〆

'の与

むわ積

値値率寄

困黒率

'''本

有有与積
園田寄黒

integerpm.pz:m <*作業用変数
cal1min'
P
ut(
X
r
)
n = irow< Xr >p = icol <Xr)
。平均値
M
E
A
l
l = ro問団 (Xr)1曲1.ω
~4 = o::~~!~x[P.p]
Md<O> = M
E
A
l
l
'
<*偏差
Xp‑
‑= x
r
‑
‑
:
‑
‑
I
::matr~x_[n ， p]*Md
C v= Xp'*Xp!dble(n?
<*分散共分散行列= di緘sq;tﾇCov.p)<*標準備差Sd= O::mãtrix[p , p]
Sd<O>= S<*変数的基準化
Xd.:Xp'CSd**(‑l))
〈事相関係数行列
R = Xd'.Xd/dble(n)EVAL= O.OdO::vector~p]
EVEC= O.OdO::matrエ玄 [p ， p]<<<<c.ll EIGEIIVV(R ， p ， EVAL~EVEC)
CUM= 1::10冒er_triJp]*EVALCRATE= EVAL/CUM[p]C聞CR =_CUM/CUMÇp~D = O.OdO::ma乞 rii: [p ， p]pm= ~基礎篇・ C 言語版

戸川隼人著・ A5 判・定価 9800 円(税込)

a~ ~0\=1.P

主要目次

ワークステーション

法の要点

C 言語の要点基礎知識線形計算

UNIX の操作

非線形方程式行列の固有値問題補間・近似・数値積分常微分方程式

icfCEVAL [i].g色 .O.999999dO) pm= i <権主成分の抽出個数i~(EVAL [il ・ E色 ~O.O ，!O)
:p~ = i <* ~より大きい固有値の数
D~i ， i] = sqrt(m世 (O.OdO ， ws))
制固有値の平方根

10
開

30
40

20

conti
FLoad = EVEC 事 D
FScore = Xd 噂 FLoad[事， 1:pm]

借主成分負荷量第( D [1 :pm ， l:pmI' 咽Ç1 :pm ， l:pm])"(・ 1】〈帯主成分得占

CBIDF=O.OdO:‑:matr x[p.
2]
do20~=O.pz-l
:
r
m
1=U
~238 j=m+1 ,pz

3 月号/発売中/定価 930 円

定叫=;，.1+1óg(~V~W)~=皿.2+EVAL [j]/ <p-m)

continue‑
0‑40 ;=1.m.

:rm.'j =rm計四2材2/(EVAL (m] -r叫〕料2

continue

曲誌唱・1・ dble(2*(p・.m)柿2句 -m+2) /dble(6*(p‑m))+r"",CHIDF[，酎 1 ， 1]=・曲*Î'm1+
nb*Cp-m) 吐。 g(rm2) <事力イ 2 乗値
CHIDF[m+1, 2] =(p‑m‑1)*(p‑m+2)/2 ‑
<*自由度

conti

write(*. 吋'棺関係数行列'

callmp:iint(R)

ite(*. 吋'固有値，寄与車，累積寄与率見方イ 2 乗値，自由度'
call 噌 rint Cl EVι; éRATE ， α皿tH. ;CHIDFI】
冒riteC+- ，*】'正規化された固有ベクトル'
call 時rint 【 EVEC)

Scheme
連載ーネットワークの辻ぱなし・証明とプログラムと・ texinfo を活用する・ Object C
でオブジェクト指向・ Mac で CLOS ・遺伝的アルゴリズムの基礎他

l
f
r
i
t
e
(
*.*】'主成分負荷量'
call mp士 int(FLo.d)

刊誌

lfrite(*. 吋， 1 以上町固有値に対応する主成分得点'

冒 rite(*: 吋以上司固有値の個数・， .
pm

ca~l :
mprintCFSc。四〉ena

低次元多様体
ゲージ理論の導入によって，低次元多様体論は大きな拡がりを経験することになりトポロジーという枠を越えて 2 ， 3 , 4 次元のゲージ理論の統ーをつくることが重要な課題になっている.
低次元多様体とはなにか/低次元物語/曲面とその写像類群/3 次元多様体の位相不変量/3次元多様体論と双曲幾何学/3 次元多様体のトポロジーと幾何学/4 次元のトポロジー他

サイエンス社
東京都千代田区神田須田町2-4 安部信ピル
宮 03-3256-1091 張替東京 7 ‑2387
1993 年 3 月号

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

(
51
)1
5
9