WANtaroHP（分析技術編）

WANtaroHP（分析技術編）　目次
01：重回帰分析（multiple regression analysis）	VB	C#	gcc
02：主成分分析（principal component analysis）	VB	C#	gcc
03：重判別分析（multiple discriminant analysis）	VB	C#	--
04：任意曲線回帰（Simplex法）	VB	C#	--
05：フーリエ級数展開	VB	C#	--

01：重回帰分析（multiple regression analysis）

概要（プログラム：vbMRA, vcsMRA, gccMRA）

重回帰分析により，回帰式：y = B0 + B1*x1 + B2*x2 +・・・+ Bm*xm の回帰係数および重相関係数を求めるプログラムです．

入力データファイル（csv形式）の並び

01 | '森北出版：多変量解析入門I_P25表1.5
02 | 4,46            
03 |  6364, 90547,19.7, 640.9,44523
04 |  7135, 12634,20.0, 128.2,10395
05 |  7266, 20101,16.3, 113.8, 8669
・・・・・

01行目：コメント
02行目：説明変数の数，データ組数
03行目：x11，x12，x13，x14，ｙ1
04行目：x21，x22，x23，x24，ｙ2
・・・以降データ組数まで繰り返し入力
注）上記は説明変数 x の数 4，データ組数を 46 とし，目的変数を y とした場合のデータ並び

出力データファイル（csv形式）の並び

01 | '森北出版：多変量解析入門I_P25表1.5
02 | y=B0+B1*x1+B2*x2+・・・Bm*xm
03 | 偏回帰係数
04 | B0=,-5508.07363049258
05 | B1=,2.31330752732194
06 | B2=,0.240910480707546
07 | B3=,-218.574370165415
08 | B4=,33.0529644259855
09 | 重相関係数：R=,0.941
10 | データ数：n=,46
11 | 残差平均：,-1.61336451444937E-11
12 | 残差標準偏差：,6357.53174194098
13 | 目的変数入力値,回帰推定値
14 | 44523,47905.2665783658
15 | 10395,13906.9412263116
16 | 8669,16341.6255537119
・・・・・

01行目：コメント
02行目：回帰式
03行目：数値の意味
04〜08行目：偏回帰係数
09〜12行目：それぞれの計算結果
13行目：数値の意味
14行目：目的変数入力値および回帰推定値
・・・以降データ組数まで繰り返し入力

使用コントロール

ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,LiatBox1

プログラムリスト

VB code	C# code	gcc code	Remarks
vb_Form1.txt	cs_Form1.txt	gccMRA.txt	プログラムリスト

入出力事例

ファイル名	概要
inpMRA5.txt	入力データ例
outMRA5.txt	出力データ例

02：主成分分析（principal component analysis）

概要（プログラム：vbPCA, vcsPCA, gccPCA）

主成分分析により固有値および固有ベクトル，主成分を求めるプログラムです．

入力データファイル（csv形式）の並び

01 | パソコンによるデータ解析入門p105
02 | 4,15
03 | 51,35,14,2
04 | 70,32,47,14
05 | 63,33,60,25
06 | 49,30,14,2
07 | 64,32,45,15
08 | 58,27,51,19
09 | 47,32,13,2
10 | 69,31,49,15
11 | 71,30,59,21
12 | 46,31,15,2
13 | 55,23,40,13
14 | 63,29,56,18
15 | 50,36,14,2
16 | 65,28,46,15
17 | 65,30,58,22

01行目：コメント
02行目：変数の数，データ組数
03行目：x11，x12，x13，x14
04行目：x21，x22，x23，x24
・・・以降データ組数まで繰り返し入力
注）上記は変数 x の数 4 ，データ組数を 15 とした場合のデータ並び

出力データファイル（csv形式）の並び

01 | パソコンによるデータ解析入門p105
02 | No.      ,1        ,2        ,3        ,4
03 | e-value  , 2.941491, 0.891700, 0.162362, 0.004448
04 | Prop     , 0.735373, 0.222925, 0.040590, 0.001112
05 | e-vector1, 0.526019, 0.319923,-0.763578, 0.194683
06 | e-vector2,-0.261656, 0.942480, 0.200582,-0.055091
07 | e-vector3, 0.578653, 0.066602, 0.227575,-0.780343
08 | e-vector4, 0.565686, 0.070322, 0.570022, 0.591717
09 | data1    ,-2.317731, 0.832810,-0.034188, 0.017762
10 | data2    , 0.900382, 0.857289,-0.662338,-0.012275
11 | data3    , 1.544318, 1.036347, 0.924547, 0.062733
12 | data4    ,-2.026758,-0.722164,-0.174767, 0.059815
13 | data5    , 0.546195, 0.639046,-0.093376, 0.007856
14 | data6    , 1.052609,-1.007354, 0.462576,-0.001419
15 | data7    ,-2.342291,-0.206211, 0.114116, 0.021857
16 | data8    , 1.051809, 0.539843,-0.545261,-0.028408
17 | data9    , 1.969426, 0.402778,-0.250594, 0.048915
18 | data10   ,-2.259044,-0.532133, 0.162490,-0.065594
19 | data11   , 0.455579,-2.392170,-0.072085, 0.026837
20 | data12   , 1.273807,-0.222464, 0.143228,-0.202216
21 | data13   ,-2.460204, 1.092547, 0.116268,-0.021841
22 | data14   , 0.966077,-0.506219,-0.421014, 0.058180
23 | data15   , 1.645825, 0.188056, 0.330397, 0.027800

01行目：コメント
02行目：番号
03行目：固有値（左から大きい順）
04行目：寄与率
05〜08行目：固有ベクトル
09〜23行目：主成分得点（標準化入力値×固有ベクトルの和）
注）１つの固有値に対応するデータは縦方向に並べて出力している

使用コントロール

ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1,ListView1

プログラムリスト

VB code	C# code	gcc code	Remarks
vb_Form1.txt	cs_Form1.txt	gccPCA.txt	プログラムリスト

入出力事例

ファイル名	概要
inpPCA.txt	分析用入力データ例
outPCA.txt	結果出力例

03：重判別分析（multiple discriminant analysis）

概要（プログラム：vbMDA, vcsMDA）

重判別分析により，入力データの正準座標，クラス別の平均正準座標を求めるプログラムです．
評価用データを入力し，正準座標と各クラス平均との Mahalanobis の汎距離を求めることもできます．
基本ルーチンは奥村晴彦先生の「パソコンによるデータ解析入門　数理とプログラミング実習」掲載のN88-BASIC(86)コードをVBに書き直したもの．
固有値を計算する Jacobi 法 Sub プロシージャは本 HP 要素技術に掲載した Jacobi 法プログラムを少し改造したものを使用しています．

入力データファイル（inpMDA.csv：csv形式）の並び

３行目以降１列目はクラス別を示します．

01 | パソコンによるデータ解析入門p105
02 | 4,15
03 | 0,51,35,14,2
04 | 0,49,30,14,2
05 | 0,47,32,13,2
06 | 0,46,31,15,2
07 | 0,50,36,14,2
08 | 1,70,32,47,14
09 | 1,64,32,45,15
10 | 1,69,31,49,15
11 | 1,55,23,40,13
12 | 1,65,28,46,15
13 | 2,63,33,60,25
14 | 2,58,27,51,19
15 | 2,71,30,59,21
16 | 2,63,29,56,18
17 | 2,65,30,58,22

01行目：コメント
02行目：変数の数，データ組数
03行目：class_1,x11，x12，x13，x14
04行目：class_2,x21，x22，x23，x24
・・・以降データ組数分まで繰り返し入力
注）上記は変数 x の数を 4，データ組数を 15 とした場合のデータ並び．03行目以降の1列目の数値は各データのクラスを示す整数で，クラス0，クラス1・・・というように0から順に入力

入出力データ例の分析用入力データは上述の通りクラスが既知でありデータ組毎にクラスを記述していますが，新規評価データ例はクラスは未知であり，判別分析により既知の各クラス平均との Mahalanobis の汎距離を算出し，入力者がクラスを判定するようにしています．
新規評価データファイルを読み込むか否かは MesageBox を表示し Yes・No を選択するようにしました．
入力データファイル（newMDA.csv：csv形式）の並びは以下の通り．クラスは未知であるため入力しません．
01 | 新規評価データ 02 | 4,6 03 | 51,35,14,2 04 | 49,30,14,2 05 | 70,32,47,14 06 | 64,32,45,15 07 | 63,33,60,25 08 | 58,27,51,19
- 01行目：コメント
- 02行目：変数の数，データ組数
- 03行目：x11，x12，x13，x14
- 04行目：x21，x22，x23，x24
- ・・・以降データ組数まで繰り返し入力
- 注）上記は変数 x の数を 4 ，データ組数を 6 とした場合のデータ並び

出力データファイル（outMDA.csv：csv形式）の並び

01 | パソコンによるデータ解析入門p105
02 | No             ,,1       ,2      ,3      ,4
03 | 固有値         ,,1513.545, 14.152,  0.000,  0.000
04 | 寄与率          ,,  0.991,  0.009,  0.000,  0.000
05 | クラス平均正準座標
06 | size=5,class0    ,-13.659,  0.377,  0.000,  0.000
07 | size=5,class1    ,  3.450, -1.333,  0.000,  0.000
08 | size=5,class2    , 10.209,  0.955,  0.000,  0.000
09 | データ毎正準座標,,1      ,2      ,3      ,4      ,dis-c0 ,dis-c1 ,dis-c2
10 | data01,class0    ,-14.596,  0.280,  0.182,  0.744,  1.214, 18.134, 24.826
11 | data02,class0    ,-13.102, -0.758,  0.688, -0.679,  1.592, 16.590, 23.394
12 | data03,class0    ,-13.836,  0.487,  0.295, -0.314,  0.478, 17.387, 24.053
13 | data04,class0    ,-12.126,  0.881, -0.874, -0.672,  1.955, 15.771, 22.362
14 | data05,class0    ,-14.638,  0.998, -0.292,  0.920,  1.508, 18.263, 24.866
15 | data06,class1    ,  2.906, -1.987, -0.433,  1.106, 16.775,  1.461,  7.962
16 | data07,class1    ,  2.612, -0.230,  0.082,  0.673, 16.297,  1.542,  7.719
17 | data08,class1    ,  4.520, -1.642, -0.600,  0.755, 18.317,  1.473,  6.327
18 | data09,class1    ,  3.233, -0.987,  0.327, -2.256, 17.100,  2.316,  7.592
19 | data10,class1    ,  3.980, -1.817,  0.624, -0.277, 17.788,  0.991,  6.852
20 | data11,class2    , 11.371,  2.866,  1.302,  0.719, 25.197,  9.087,  2.686
21 | data12,class2    ,  8.361,  1.142,  0.200, -1.113, 22.063,  5.615,  2.175
22 | data13,class2    , 10.309, -1.106,  0.601,  0.565, 24.028,  6.912,  2.223
23 | data14,class2    , 10.340,  1.032, -2.664, -0.290, 24.157,  7.761,  2.684
24 | data15,class2    , 10.665,  0.842,  0.561,  0.119, 24.336,  7.557,  0.741
25 | 新規評価データ正準座標
26 | new   ,00        ,-14.596,  0.280,  0.182,  0.744,  1.214, 18.134, 24.826
27 | new   ,01        ,-13.102, -0.758,  0.688, -0.679,  1.592, 16.590, 23.394
28 | new   ,02        ,  2.906, -1.987, -0.433,  1.106, 16.775,  1.461,  7.962
29 | new   ,03        ,  2.612, -0.230,  0.082,  0.673, 16.297,  1.542,  7.719
30 | new   ,04        , 11.371,  2.866,  1.302,  0.719, 25.197,  9.087,  2.686
31 | new   ,05        ,  8.361,  1.142,  0.200, -1.113, 22.063,  5.615,  2.175

01行目：コメント
02行目：番号
03行目：固有値（左から大きい順）
04行目：寄与率
05行目：数値の意味
06〜08行目：データ組数（size），クラス別とそれぞれの固有値に対応する平均正準座標
09行目：数値の意味
10〜24行目：データ毎のクラス別と正準座標．dis-c0〜dis-c2 はクラス 0〜2 の平均正準座標からの距離
25行目：数値の意味
26〜31行目：新規評価データの正準座標．7〜9列は dis-c0〜dis-c2 の値であり，クラス0〜2 の平均正準座標からの距離
注）新規評価データについては dis-c0〜dis-c2 の値を見て最も近い近い距離にあるクラスに属するものと評価すればよい

使用コントロール

ToolStrip1(Button1),OpenFileDialog1,SaveFileDialog1

プログラムリスト

VB code	C# code	gcc code	Remarks
vb_Form1.txt	cs_Form1.txt	----	プログラムリスト

入出力事例

ファイル名	概要
inpMDA.txt	分析用入力データ例
newMDA.txt	新規評価データ例
outMDA.txt	結果出力例

04：任意曲線回帰（Simplex法）

概要（プログラム：vbSIMPLEX, vcsSIMPLEX）

シンプレックス法により任意曲線への回帰式を求めるプログラムです．
基本ルーチンは奥村晴彦先生の「パソコンによるデータ解析入門　数理とプログラミング実習」掲載の N88-BASIC(86) コードを VB および C# に書き直したものであり，データプロットと求めた回帰式の関係を確認するためグラフ出力ルーチンを加えています．
このプログラムでは，データファイルの他，プログラム本体についても，回帰曲線が変わるたびに，回帰係数の数（パラメータ数，変数：MPARA）を指定する Function プロシージャ（MPARA_SET），回帰式を定義する Function プロシージャ（VFunc）を書き直す必要があります．
また回帰係数の初期値・検索ステップ・収束判定基準を定義する Sub プロシージャ INIVAL も必要に応じて書き直した方が効率がよいと思われます（逆にまとはずれの数値だと収束しない可能性がある）．
しかしながら初期値の工夫など試行錯誤はある程度必要ですが，回帰できない確率は低く重宝するプログラムです．
このプログラムでは Form1 に回帰結果のグラフを示し，MessageBox に回帰係数を示すようにしました．
実行すると，回帰係数を書き出す csv ファイルと回帰結果のグラフを書き出す png 画像ファイルが出力されます．

入力データファイル（csv形式）のデータ並び

01 | パソコンによるデータ解析入門p236,
02 | エネルギー,反応率
03 | 600,70,20
04 | 400,50,20
05 | -10.0,10.0,2.0
06 | 0.0,10.0,2.0
07 | 11
08 | -10,2.127
09 | -8,2.52
10 | -6,2.629
11 | -4,2.938
12 | -2,3.414
13 | 0,4.669
14 | 2,8.014
15 | 4,6.372
16 | 6,4.596
17 | 8,4.296
18 | 10,4.291

01行目：コメント
02行目：x 軸名（単位含む），y 軸名（単位含む）
03行目：画像横サイズ（ピクセル），画像左余白（縦軸描画スペース：ピクセル），画像右余白（ピクセル）
04行目：画像縦サイズ（ピクセル），画像下余白（横軸描画スペース：ピクセル），画像上余白（ピクセル）
05行目：x 軸最小値，x 軸最大値，x 軸増分
06行目：y 軸最小値，y 軸最大値，y 軸増分
07行目：データ組数
08行目以降：x 値，y 値・・・・以降同様
注）２行目から６行目までは２次元グラフ画像出力プログラムと同じ考え方の入力．ただしここでのグラフ出力は普通軸のみを対象としています．

出力データファイル（csv形式）のデータ並び

01 | パソコンによるデータ解析入門p236
02 | iflag=,0
03 | iteration=,412
04 | MPARA=,5
05 | parameter(0)=,2.522624E+000
06 | parameter(1)=,1.656142E+000
07 | parameter(2)=,5.242037E+000
08 | parameter(3)=,3.054768E+000
09 | parameter(4)=,9.810139E-002

01行目：コメント
02行目：iflag=0 なら収束，iflag=1 なら未収束
03行目：iterasion は収束までの繰り返し回数を示す
04行目：MPARA は回帰係数の数を示す
05〜09行目：回帰係数の値

入出力データ事例について

入力データ例としては，例１として前出奥村先生の著書に掲載されているローレンツ型反応曲線＋１次関数への回帰，例２としてバラツキを加味して作成した減衰振動曲線への回帰を示しています．

使用コントロール

ToolStrip1(Button1),PictureBox1,PictureBox2,OpenFileDialog1,SaveFileDialog1

プログラムリスト

VB code	C# code	gcc code	Remarks
vb_Form1.txt	cs_Form1.txt	----	プログラムリスト

入出力事例

ファイル名	概要
inpSIMP0.txt	分析用入力データ例（例１）
outSIMP0.txt	結果出力例（例１）
outSIMP0.png	結果出力図例（例１）
inpSIMP1.txt	分析用入力データ例（例２）
outSIMP1.txt	結果出力例（例２）
outSIMP1.png	結果出力図例（例２）
pdfSIMP.pdf	例１および例２のグラフ出力

05：フーリエ級数展開

概要（プログラム：vbFAPP, vcsFAPP）

フーリエ級数展開により任意曲線への回帰推定値を求めるプログラムです．
回帰式：y=f(x)=A0/2+Σ[Ak*cos(2πk*x/N/Δx)+Bk*sin(2πk*x/N/Δx)] k=1,・・・, mm
ここに N はデータ数（偶数個），Δx は x 値の刻み（一定値）です．π はパイ（円周率）です．
級数和をとる項数 mm は本来，k=1 から k=N/2-1 までですが，ここでは mm までとして近似しています．
フーリエ係数 Ak および Bkは，高速フーリエ変換により求めています．Ak=2*Re(Ck), Bk=-2*Im(Ck) で，Ck は複素フーリエ係数です．すなわち Ak は複素フーリエ係数の実数部の 2 倍，Bk は複素フーリエ係数の虚数部の 2 倍値にマイナス符号をつけたものです．
高速フーリエ変換ではデータ数を 2 の累乗にするため「後続の 0」を入力データの末尾に足しています．このため入力データの最後の y 値が 0 となるように y 軸を平行移動させ，フーリエ係数を求めた後に再度 y 軸を平行移動させて元の入力データとの整合をとっています．
x 軸のデータは等間隔で入力する必要があります．
級数和をとる項数 mm はデータファイルで入力します．
実行すると，回帰結果を書き出す csv ファイルと回帰結果のグラフを書き出す png 画像ファイルが出力されます．
図中の赤丸は入力データ，青線が回帰推定値，緑丸は残差です．
ここでの入出力事例では，級数和をとる項数は，mm=5 と mm=100 としています．

入力データファイル（csv形式）のデータ並び

01 | sample_mm=5
02 | 横軸ｘ,縦軸ｙ
03 | 1000,70,20
04 | 500,50,30
05 | 1000,8000,1000
06 | -5,30,5
07 | 6299,5
08 | 1000,11.3
09 | 1001,11.2
10 | 1002,11.1
11 | 1003,11.0
・・・・

01行目：コメント
02行目：x 軸名（単位含む），y 軸名（単位含む）
03行目：画像横サイズ（ピクセル），画像左余白（縦軸描画スペース：ピクセル），画像右余白（ピクセル）
04行目：画像縦サイズ（ピクセル），画像下余白（横軸描画スペース：ピクセル），画像上余白（ピクセル）
05行目：x 軸最小値，x 軸最大値，x 軸増分
06行目：y 軸最小値，y 軸最大値，y 軸増分
07行目：データ組数，級数和をとる項数
08行目以降：x 値，y 値・・・・以降同様
注）２行目から６行目までは２次元グラフ画像出力プログラムと同じ考え方の入力．ただしここでのグラフ出力は普通軸のみを対象としています．

出力データファイル（csv形式）のデータ並び

01 | sample_mm=5
02 | ndata,6299
03 | x,y,func
04 | 1.000000E+003,1.130000E+001,1.475526E+001
05 | 1.001000E+003,1.120000E+001,1.475746E+001
06 | 1.002000E+003,1.110000E+001,1.475967E+001
07 | 1.003000E+003,1.100000E+001,1.476188E+001
・・・・

01行目：コメント
02行目：データ組数
03行目：出力数値並びの意味
04行目以降：x 値，y 値，回帰推定値

使用コントロール

ToolStrip1(Button1),PictureBox1,PictureBox2,OpenFileDialog1,SaveFileDialog1

プログラムリスト

VB code	C# code	gcc code	Remarks
vb_Form1.txt	cs_Form1.txt	----	プログラムリスト

入出力事例

ファイル名	概要
inp_mm005.txt	分析用入力データ例（例１）
out_mm005.txt	結果出力例（例１）
fig_out_mm005.png	結果出力図例（例１）
inp_mm100.txt	分析用入力データ例（例２）
out_mm100.txt	結果出力例（例２）
fig_out_mm100.png	結果出力図例（例２）