笛田薫Blog: 12月 2006

2006-12-30

汐留シオサイト

汐留シオサイトの日本テレビビルに、宮崎駿さんデザインの大時計がありました。12月21日に公開されたばかりだそうです。

2006-12-29

アクアミュージアム

横浜八景島シーパラダイスアクアミュージアムへ行ってきました。

綺麗な魚が出迎えしてくれました。

ペンギンは水の中を飛ぶように泳ぎます。

海亀。

イルカの水槽の下を通れます。

ペリカンに餌をあげることが出来ます。

口が大きく広がりますね。

2006-12-26

2地点の緯度、経度から距離を求める

球面三角法とその角度
これらを使って2地点の緯度、経度から間の角度を求めることが出来ます。
但し、地球は完全な球ではなく楕円球なので補正が必要になります。

間の角度をラジアンで求めたら、それに地球の半径をかけて距離を求めることが出来ます。

手っ取り早く求めるには
http://www2.neweb.ne.jp/wd/nobuaki/New_Homepage/okinawa703.htm
http://oshiete.eibi.co.jp/kotaeru.php3?q=249931
http://www1.ths.titech.ac.jp/club/sci_club/astronomy/ISSP/Cal/HENKAN.html
あたりが参考になりそう。

2006-12-21

Excel VBAの解答例

今日はパワーポイントを勉強しますが、その前に前回のレポート課題である、Visual Basicを使った平均、分散を求めるプログラムの解答例を示します。

Sub 解答例()
Dim x()

i = 1
Do While Cells(i, 1) <> ""
    i = i + 1
Loop
kosu = i - 1
Cells(1, 2) = "個数"
Cells(1, 3) = kosu
ReDim x(kosu)

For i = 1 To kosu
    x(i) = Cells(i, 1)
Next i

wa = 0
For i = 1 To kosu
    wa = wa + x(i)
Next i
Cells(2, 2) = "合計"
Cells(2, 3) = wa

heikin = wa / kosu
Cells(3, 2) = "平均"
Cells(3, 3) = heikin

wa = 0
For i = 1 To kosu
    wa = wa + (x(i) - heikin) ^ 2
Next i
Cells(4, 2) = "分散"
Cells(4, 3) = wa / kosu

End Sub

2006-12-20

第13回講義:連立一次方程式第3回

式の数、未知数の数が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

のように書けないでしょうか。そうすれば繰り返しの回数を増やすことで式の数が増えても解けるようになります。

2006-12-18

環境データ解析の例

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

都市番号	都市	気温y	緯度x₁	経度x₂	標高x₃
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₃とおきます。緯度、経度、標高と気温の関係は色々考えることが出来ますが、今回は緯度、経度、標高の3変量を使う最も簡単なものとして

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以上です。だからグラフは下に凸になっていて、どこかで最小値を取ります。

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

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

β₀

β₁

β₂

β₃

をΒとおき、

ベクトル

y₁

y₂

・

y₂₀

をYとおき、

最後に行列

1	x₁₁	x₁₂	x₁₃
1	x₂₁	x₂₂	x₂₃
・	・	・	・
・	・	・	・
1	x_20,1	x_20,2	x_20,3

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

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

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

この連立一次方程式を解いて、β₀,β₁,β₂,β₃を求めてください。正解は、 β₀=38.3、 β₁=-1.17、 β₂=0.02、 β₃=-0.0098、そして神戸の気温は、38.3-1.17×34.68+0.02×135.18-0.0098×59.3=-0.15です。実際の気温は1.2でした。

2006-12-13

第12回講義:連立一次方程式第2回

このテーマは奥が深いので、今後の課題を紹介しておきます。まずは焦らず前回の課題に取り組んでください。

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

今回はその準備として、似たような処理を繰り返し書いている部分をFor文を使って1回書くだけで済ませましょう。

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 (2,1)成分を0にするために第二行から第一行の(2,1)成分倍を引くプログラムを書いてください。
Rem (3,1)成分を0にするために第三行から第一行の(3,1)成分倍を引くプログラムを書いてください。

Rem (2,2)成分を1にするために第二行を(2,2)成分で割るプログラムを書いてください。

Rem (1,2)成分を0にするために第一行から第二行の(1,2)成分倍を引くプログラムを書いてください。
Rem (3,2)成分を0にするために第三行から第二行の(3,2)成分倍を引くプログラムを書いてください。

Rem (3,3)成分を1にするために第三行を(3,3)成分で割るプログラムを書いてください。

Rem (1,3)成分を0にするために第一行から第三行の(1,3)成分倍を引くプログラムを書いてください。
Rem (2,3)成分を0にするために第二行から第三行の(2,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箇所の色をつけた部分は、それぞれ2回ずつ殆ど似たようなことを書いていますよね。これをFor文でまとめたいと思います。

Rem (2,1)成分を0にするために第二行から第一行の(2,1)成分倍を引くプログラムを書いてください。
Rem (3,1)成分を0にするために第三行から第一行の(3,1)成分倍を引くプログラムを書いてください。

は

For p=2 to 3
Rem (p,1)成分を0にするために第p行から第一行の(p,1)成分倍を引くプログラムを書いてください。
Next p

と書けば良いですし、

Rem (1,3)成分を0にするために第一行から第三行の(1,3)成分倍を引くプログラムを書いてください。
Rem (2,3)成分を0にするために第二行から第三行の(2,3)成分倍を引くプログラムを書いてください。

は

For p=1 to 2
Rem (p,3)成分を0にするために第p行から第三行の(p,3)成分倍を引くプログラムを書いてください。
Next p

と書けば良いです。では

Rem (1,2)成分を0にするために第一行から第二行の(1,2)成分倍を引くプログラムを書いてください。
Rem (3,2)成分を0にするために第三行から第二行の(3,2)成分倍を引くプログラムを書いてください。

の部分はどうすればよいか、考えてみてください。

2006-12-10

故郷

3年ぶりに実家に帰ってきました。

久しぶりに帰る家は、一歩入った瞬間に「あれ、こんなに天井低かったっけ?電気のスイッチはこんなに低かったかな?」と戸惑いました。3年前と比べて背が伸びているわけではないですが(笑)、記憶の中の家は3年前ではなく、子供の頃の視線からの記憶だからでしょう。
遠くに感じた小学校も、懐かしさを感じながら歩くと近く感じますし、ちょっとした探検気分で登ったその向こうの山も、切り開かれて道路が通っていました。
逞しかった親父の足も細くなって(それは僕がかじりすぎたから?)
次に帰るのはいつになるんだろう。

福岡も大きく変わっていました。バブル崩壊で痛手を受けた商業ビルも、新しい経営者に変わって生まれ変わっていました。
また帰りたいなぁ。

2006-12-07

情報基盤センターでのVBA

情報基盤センターでVBAを使うときのコツが掴めてきました。
VBAがくっついたExcelファイルを一旦保存して、改めて開くとプログラムが動かないのですから、動かないファイルを開いた状態でもう一つExcelのbookを新規作成して、そちらにプログラムもExcelシートのデータもコピーすれば動きます。

このブログで表を書く方法もコツが掴めてきました。どうやらここのHTMLは、文中の改行が無視されないらしく、TABLEのタグの途中で改行すると、改行した個数だけ表の前に空行が空きます。
正確に言うと、＜ＴＤ＞と＜／ＴＤ＞の間、つまり表の中に書かれる部分で改行すると表示される表もその部分で改行され、それ以外の部分で改行すると表の前で改行されます。
従って、表の中に空の行がある場合には、その行の＜ＴＤ＞と＜／ＴＤ＞の間で一回だけ改行しておかないと行が潰れますし、逆に改行して欲しくないけれどHTMLのソースを見やすくするためにソースでは改行したい場合は＜！－－ここで改行
－－＞と書いて、改行をコメントアウトしてしまう必要があります。

さらに、このことを書いていて分かりましたが、HTMLでは無い方の、WYSIWYG画面でもタグと同じ文字を書くとタグと認識されてしまいます。よって上記説明は全角文字を使いました。

横幅も広く取れるようになりましたし、これで講義の説明用に使えそうです。

2006-12-06

第11回講義:連立一次方程式第1回

前回の解答

Sub 掛算()
Dim x() As Single, y() As Single, z() As Single
Dim gyou1 As Integer, gyou2 As Integer, retsu1 As Integer, retsu2 As Integer
Dim 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
retsu1 = 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
retsu2 = j - 1

If retsu1 = gyou2 Then
　ReDim x(gyou1, retsu1)
　ReDim y(gyou2, retsu2)
　ReDim z(gyou1, retsu2)

　Sheets(1).Select
　For i = 1 To gyou1
　　For j = 1 To retsu1
　　　x(i, j) = Selection.Cells(i, j)
　　Next j
　Next i

　Sheets(2).Select
　For i = 1 To gyou2
　　For j = 1 To retsu2
　　　y(i, j) = Selection.Cells(i, j)
　　Next j
　Next i

　For i = 1 To gyou1
　　For j = 1 To retsu2
　　　z(i, j) = 0
　　　For k = 1 To retsu1
　　　　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 retsu2
　　　Selection.Cells(i, j) = z(i, j)
　　Next j
　Next i

Else
　Sheets(3).Select
　Selection.Cells(1, 1) = "ベクトルの型が異なるので掛け算できません"
End If

End Sub

履き出し法

この問題をExcelを使って解きましょう。まず係数の部分と定数項の部分を Excelのセルに入力します。

A	B	C
2	4	4
2	-1	-6
2	3	-2
Sheet1

A	B	C
18
-2
-6
Sheet2

まずこの係数部分と定数項を配列に読み込んで、Sheets(3)の(1,1)セルから順に表示します。

Sub 連立一次方程式1()
Dim x(3, 4) As Single
Dim i As Integer, j 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
Cells(1, 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

するとこのようになります。

A	B	C	D
2	4	4	18
2	-1	-6	-2
2	3	-2	-6
Sheet3

履き出し法の実行自分で計算する手順を考えながら、プログラムを書きます。途中経過もExcelに出力すると分かりやすいです。

A	B	C	D
2	4	4	18
2	-1	-6	-2
2	3	-2	-6

1	2	2	9	(1,1)成分を1にするために第1行を(1,1)成分で割ります
2	-1	-6	-2
2	3	-2	-6

1	2	2	9
0	-5	-10	-20	(2,1)成分を0にするために第2行から第1行×(2,1)成分を引きます
0	-1	-6	-24	(3,1)成分を0にするために第3行から第1行×(3,1)成分を引きます
Sheet3

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)成分で割るプログラムを書いてください。

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

Rem (2,1)成分を0にするために第二行から第一行の(2,1)成分倍を引くプログラムを書いてください。
Rem (3,1)成分を0にするために第三行から第一行の(3,1)成分倍を引くプログラムを書いてください。

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

演習問題

1. 上のプログラムにおいて、Rem 文で書いている手順をBasicのプログラムで書きなさい。
最後まで解くとこうなります。

A	B	C	D	E(ここは説明です)
2	4	4	18
2	-1	-6	-2
2	3	-2	-6

1	2	2	9	(1,1)成分を1にするために第1行を(1,1)成分で割ります
2	-1	-6	-2
2	3	-2	-6

1	2	2	9
0	-5	-10	-20	(2,1)成分を0にするために第2行から第1行×(2,1)成分を引きます
0	-1	-6	-24	(3,1)成分を0にするために第3行から第1行×(3,1)成分を引きます

1	2	2	9
0	1	2	4	(2,2)成分を1にするために第2行を(2,2)成分で割ります
0	-1	-6	-24

1	0	-2	1	(1,2)成分を0にするために第1行から第2行×(1,2)成分を引きます
0	1	2	4
0	0	-4	-20	(3,2)成分を0にするために第3行から第2行×(3,2)成分を引きます

1	0	-2	1
0	1	2	4
0	0	1	5	(3,3)成分を1にするために第3行を(3,3)成分で割ります

1	0	0	11	(1,3)成分を0にするために第1行から第3行×(1,3)成分を引きます
0	1	0	-6	(2,3)成分を0にするために第2行から第3行×(2,3)成分を引きます
0	0	1	5
Sheet3

2006-12-05

商圏分析

岡山県の平成16年度商業統計調査確報
第３表　小売業　市町村別　売場面積規模別事業所数・従事者数・年間商品販売額・売場面積
ここにExcel形式で置いてある。
データは
(財)統計情報研究開発センター(Sinfonica)から入手。

昨日は来年度のゼミ生が、現在のゼミを見学。
希望は地理データ、数理統計。一人は就活で欠席。

夜は並列計算の使い方。
Cross Validationやbootstrapのように、互いに依存関係の無い計算を多数回繰り返すときに、それぞれの計算機に分けて計算させて、最後に集計させるようにすれば良い。

2006-12-04

連立方程式の解法

連立方程式の解法の最終回。
ピボットが0になるときにどう対処するか、ということで前回質問がありました。
特異値分解からMoore-Penrose一般逆行列を求める、という方向を詳説していると講義の本筋から逸れますので、この講義では変数選択で対応することにして、
ピボットが極めて0に近い→答の精度が低くてもとりあえず求める。
ピボットが本当に0でOverflowになる→0で割った商の代わりに0そのものを用いる。
という方法で対処することにしました。

一般逆行列に興味がある人は、例えば東大出版会の線型代数演習(齋藤正彦著)を見てください。

2006-12-03

インターネット・マーケティング

来年のゼミ生の中に、マーケティングに興味があるという学生さんが居たので、参考になるかも。
Googleアドセンス情報日本語非公式版のインターネット・マーケティングのおすすめ資料
THE Google Adsense

2006-12-02

論文の戻し

某誌に投稿していた某論文がone “Accept after minor revision” and two “Reconsider after major revision”で戻ってきました。reviseせねば。

2006-12-01

ポアソン分布

重要な分布として、二項分布やポアソン分布を教えました。
ポアソン分布は二項分布の近似として教えたのですが、教科書の例「大学の全学生のうちで1年間に死亡する学生の数Xはポアソン分布に従うものと考えられる。」は、もう少し正確に言うと「二項分布に従うけれど、nが大きくpの値が極めて小さいのでポアソン分布で近似できる。」と言うべきでしょうね。
あまり堅い話になると、本当に二項分布なのかも怪しいです。突発事故で死亡者が多数出ることを考えると独立が怪しいので。

話は変わりますが、入門用に分かりやすく説明したページを見つけました。
早稲田大学人間科学部　向後千春研究室

登録: 投稿 (Atom)

笛田薫Blog

2006-12-30

汐留シオサイト

2006-12-29

アクアミュージアム

2006-12-26

2地点の緯度、経度から距離を求める

2006-12-21

Excel VBAの解答例

2006-12-20

第13回講義:連立一次方程式第3回

2006-12-18

環境データ解析の例

2006-12-13

第12回講義:連立一次方程式第2回

2006-12-10

故郷

2006-12-07

情報基盤センターでのVBA

2006-12-06

第11回講義:連立一次方程式第1回

2006-12-05

商圏分析

2006-12-04

連立方程式の解法

2006-12-03

インターネット・マーケティング

2006-12-02

論文の戻し

2006-12-01

ポアソン分布

Translate

最近１カ月の人気の投稿

ラベルの説明

アクセス数

マイブログリスト

過去の日記

最近のコメント

2006-12-30

2006-12-29

2006-12-26

2006-12-21

2006-12-20

2006-12-18

2006-12-13

2006-12-10

2006-12-07

2006-12-06

2006-12-05

2006-12-04

2006-12-03

2006-12-02

2006-12-01

Translate

最近１カ月の人気の投稿

ラベルの説明

アクセス数

マイブログ リスト

過去の日記

最近のコメント

マイブログリスト