連立一次方程式の解法

1. 前回の課題の解答

   まず2行2列の行列の積のプログラムです

   Sub 行列の掛算1()
   Dim a(2, 2), b(2, 2), c(2, 2)
   
   Sheets(1).Select
   For i = 1 To 2
   　 For j = 1 To 2
   　　　 a(i, j) = Selection.Cells(i, j)
   　 Next j
   Next i
   
   Sheets(2).Select
   For i = 1 To 2
   　 For j = 1 To 2
   　　　 b(i, j) = Selection.Cells(i, j)
   　 Next j
   Next i
   
   For i = 1 To 2
   　 For j = 1 To 2
   　　　 c(i, j) = 0
   　　　 For k = 1 To 2
   　　　　　 c(i, j) = c(i, j) + a(i, k) * b(k, j)
   　　　 Next k
   　　Next j
   Next i
   
   Sheets(3).Select
   For i = 1 To 2
   　 For j = 1 To 2
   　　　 Selection.Cells(i, j) = c(i, j)
   　 Next j
   Next i
   End Sub

   次に行列の型を調べて掛け算できるかどうか判断するプログラムです。

   Sub 行列の掛算2()
   Dim x() As Single, y() As Single, z() As Single, gyou1 As Integer, gyou2 As Integer, retu1 As Integer, retu2 As Integer, i As Integer, j As Integer, k As Integer
   
   Sheets(1).Select
   i = 1
   Do While Selection.Cells(i, 1) <> "" And IsNumeric(Selection.Cells(i, 1))
   　i = i + 1
   Loop
   gyou1 = i - 1
   
   j = 1
   Do While Selection.Cells(1, j) <> "" And IsNumeric(Selection.Cells(1, j))
   　j = j + 1
   Loop
   retu1 = j - 1
   
   Sheets(2).Select
   i = 1
   Do While Selection.Cells(i, 1) <> "" And IsNumeric(Selection.Cells(i, 1))
   　i = i + 1
   Loop
   gyou2 = i - 1
   
   j = 1
   Do While Selection.Cells(1, j) <> "" And IsNumeric(Selection.Cells(1, j))
   　j = j + 1
   Loop
   retu2 = j - 1
   
   If retu1 = gyou2 Then
   　ReDim x(gyou1, retu1)
   　ReDim y(gyou2, retu2)
   　ReDim z(gyou1, retu2)
   　
   　Sheets(1).Select
   　For i = 1 To gyou1
   　　For j = 1 To retu1
   　　　x(i, j) = Selection.Cells(i, j)
   　　Next j
   　Next i
   　
   　Sheets(2).Select
   　For i = 1 To gyou2
   　　For j = 1 To retu2
   　　　y(i, j) = Selection.Cells(i, j)
   　　Next j
   　Next i
   　
   　For i = 1 To gyou1
   　　For j = 1 To retu2
   　　　z(i, j) = 0
   　　　For k = 1 To retu1
   　　　　z(i, j) = z(i, j) + x(i, k) * y(k, j)
   　　　Next k
   　　Next j
   　Next i
   　
   　Sheets(3).Select
   　For i = 1 To gyou1
   　　For j = 1 To retu2
   　　　Selection.Cells(i, j) = z(i, j)
   　　Next j
   　Next i
   　
   Else
   　Msgbox ("行列の型が合わないので掛け算できません")
   End If
   
   End Sub

2. 一般化

   式の数、未知数の数が4個になると前回のプログラムは使えません。 折角作ったプログラムですので、式の数が増えても使えるようにしましょう。

   Sub 連立一次方程式2() Dim x(3, 4) As Single Dim i As Integer, j As Integer, basho As Integer Sheets(1).Select For i = 1 To 3 For j = 1 To 3 x(i, j) = Selection.Cells(i, j) Next j Next i Sheets(2).Select For i = 1 To 3 x(i, 4) = Selection.Cells(i, 1) Next i Sheets(3).Select basho = 1 Cells(basho, 1).Select For i = 1 To 3 For j = 1 To 4 Selection.Cells(i, j) = x(i, j) Next j Next i Rem (1,1)成分を1にするために第一行を(1,1)成分で割るプログラムを書いてください。 Rem p=2,3に対して、(p,1)成分を0にするために第p行から第一行の(p,1)成分倍を引くプログラムを書いてください。 Rem (2,2)成分を1にするために第二行を(2,2)成分で割るプログラムを書いてください。 Rem p=1,3に対して、(p,2)成分を0にするために第p行から第二行の(p,2)成分倍を引くプログラムを書いてください。 Rem (3,3)成分を1にするために第三行を(3,3)成分で割るプログラムを書いてください。 Rem p=1,2に対して、(p,3)成分を0にするために第p行から第三行の(p,3)成分倍を引くプログラムを書いてください。 basho = basho + 4 Cells(basho, 1).Select For i = 1 To 3 For j = 1 To 4 Selection.Cells(i, j) = x(i, j) Next j Next i End Sub

   このプログラムの赤い部分、青い部分、緑の部分は非常に似ています。このように似た部分を3回書くのではなく
   For q=1 to 3
   　　共通の部分
   Next q
   のように書けないでしょうか。そうすれば繰り返しの回数を増やすことで式の数が増えても解けるようになります。

3. データ解析

   未知数の数が4つの連立一次方程式も解ける様になった人は次の課題に挑戦してください。

   以下のように日本の20都市に対して、１月の日最低気温の月平均値、緯度、経度、標高のデータがあります。
   

   都市番号	都市	気温y	緯度x1	経度x2	標高x3
   1	稚内	-8.0	45.42	141.68	2.8
   2	旭川	-13.6	43.77	142.37	111.9
   3	札幌	-9.5	43.05	141.33	17.2
   4	青森	-5.4	40.82	140.78	3.0
   5	盛岡	-6.7	39.70	141.17	155.2
   6	仙台	-3.2	38.27	140.90	38.9
   7	金沢	-0.1	36.55	136.65	26.1
   8	長野	-5.5	36.67	138.20	418.2
   9	高山	-7.6	36.15	137.25	560.2
   10	軽井沢	-10.0	36.33	138.55	999.1
   11	名古屋	-0.9	35.17	136.97	51.1
   12	飯田	-4.7	35.52	137.83	481.8
   13	東京	-0.4	35.68	139.77	5.3
   14	鳥取	0.5	35.48	134.23	7.1
   15	京都	-0.6	35.02	135.73	41.4
   16	広島	0.2	34.37	132.43	29.3
   17	福岡	1.5	33.58	130.38	2.5
   18	鹿児島	2.0	31.57	130.55	4.3
   19	高知	0.1	33.55	133.53	1.9
   20	那覇	13.5	26.23	127.68	34.9

   出典:坂元慶行、石黒真木夫、北川源四郎著「情報量統計学」共立出版株式会社
   
   このデータを元に、緯度、経度、標高と気温の関係を数式で表したいと思います。 また神戸の緯度は34.68、経度は135.18、標高は59.30です。神戸の気温を推定しましょう。

   気温をy, 緯度をx₁、経度をx₂、標高をx₃とおきます。 緯度、経度、標高と気温の関係は色々考えることが出来ますが、ここでは簡単なものとして

   y=β₀+β₁x₁+β₂x₂+β₃x₃+誤差

   という形を考えます。
   β₀,β₁,β₂,β₃はどうやって決めましょうか。
   本当の気温yと、緯度、経度、標高から推定した気温β₀+β₁x₁+β₂x₂+β₃x₃ との差が小さくなるように決めましょう。

   データの
   1番目の都市の気温をy₁、緯度をx₁₁、経度をx₁₂、標高をx₁₃
   2番目の都市の気温をy₂、緯度をx₂₁、経度をx₂₂、標高をx₂₃
   ・・・・・
   20番目の都市の気温をy₂₀、緯度をx_20,1、経度をx_20,2、標高をx_20,3
   と書きます。すると、本当の気温と推定した気温の差は

   y₁-(β₀+β₁x₁₁+β₂x₁₂+β₃x₁₃)
   y₂-(β₀+β₁x₂₁+β₂x₂₂+β₃x₂₃)
   ・・・・・
   y₂₀-(β₀+β₁x_20,1+β₂x_20,2+β₃x_20,3)

   です。これらはプラスになったりマイナスになったりするので、そのまま足すと打ち消しますから、二乗して足して

   ∑_i=1²⁰　{y_i-(β₀+β₁x_i1+β₂x_i2+β₃x_i3)}²

   が最小になるようなβ₀,β₁,β₂,β₃を求めることにします。
   一見、求めるのは難しそうに見えます。でもこれはβ₀,β₁,β₂,β₃に関する二次式で、しかもは二乗して足していますから0以上です。だからグラフは下に凸になっていて、どこかで最小値を取ります。

   公式の導出は省略しますが、実は最小にするβ₀,β₁,β₂,β₃は次の連立一次方程式を解くことで求めることが出来ます。

   まず未知数β₀,β₁,β₂,β₃を縦に並べた ベクトル

   β0

   β1

   β2

   β3

   をΒとおき、
   
   ベクトル

   y1

   y2

   ・

   ・

   y20

   をYとおき、
   
   最後に行列

   1	x11	x12	x13
   1	x21	x22	x23
   ・	・	・	・
   ・	・	・	・
   1	x20,1	x20,2	x20,3

   をXとおきます。
   1が縦に並んでいるのはβ₀の係数だからです。この時、最小にするβ₀,β₁,β₂,β₃は 連立一次方程式

   (X^t・X)Β=(X^t・Y)

   を解くことで得られます。Xの右上のtは転置を表します。行列X、ベクトルYの配置は上のデータの表と同じですからコピーすれば入力しやすいです。

   この連立一次方程式を解いて、β₀,β₁,β₂,β₃を求めてください。