では、マラリア患者を刺した蚊の中で、人に感染させられる蚊の割合は？
話を簡単にするために、人に感染させられるようになるまでに要する日数を10日に固定すると…
日生存率をrとすると、刺して0日経過から9日経過の数は、刺した蚊の数で割って基準化すると
1+r+r2+…+r9
一方感染力を持つ蚊はr10+r11+…
よって比率は
r10+r11+…
－－－－－－－－－－－－－－－－－=r10
1+r+r2+…+r9+r10+r11+…
なんか、凄く簡単なことを長々計算したような…
| index | city | 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 |
緯度、経度、標高と気温の関係は色々考えることが出来ますが、今回は緯度、経度、標高の3変量を使う最も簡単なものとして
y=β0+β1x1+β2x2+β3x3+誤差
という形を考えます。
β0,β1,β2,β3はどうやって決めましょうか。
本当の気温yと、緯度、経度、標高から推定した気温β0+β1x1+β2x2+β3x3
との差が小さくなるように決めましょう。
データの
1番目の都市の気温をy1、緯度をx11、経度をx12、標高をx13
2番目の都市の気温をy2、緯度をx21、経度をx22、標高をx23
・・・・・
20番目の都市の気温をy20、緯度をx20,1、経度をx20,2、標高をx20,3
と書きます。すると、本当の気温と推定した気温の差は
y1-(β0+β1x11+β2x12+β3x13)
y2-(β0+β1x21+β2x22+β3x23)
・・・・・
y20-(β0+β1x20,1+β2x20,2+β3x20,3)
です。これらはプラスになったりマイナスになったりするので、そのまま足すと打ち消しますから、二乗して足して
∑i=120 {yi-(β0+β1xi1+β2xi2+β3xi3)}2
が最小になるようなβ0,β1,β2,β3を求めることにします。
一見、求めるのは難しそうに見えます。でもこれはβ0,β1,β2,β3に関する二次式で、しかもは二乗して足していますから0以上です。だからグラフは下に凸になっていて、どこかで最小値を取ります。
公式の導出は書きにくいので講義のときに板書しますが、実は最小にするβ0,β1,β2,β3は次の連立一次方程式を解くことで求めることが出来ます。
まず未知数β0,β1,β2,β3を縦に並べたベクトル
| β0 |
| β1 |
| β2 |
| β3 |
| y1 |
| y2 |
| ・ |
| ・ |
| y20 |
| 1 | x11 | x12 | x13 |
| 1 | x21 | x22 | x23 |
| ・ | ・ | ・ | ・ |
| ・ | ・ | ・ | ・ |
| 1 | x20,1 | x20,2 | x20,3 |
ハフモデルについて詳しくはRetail Trade Area Analysis Using the Huff Model - Articlesをご覧ください。
まず正しいコインを10回投げたときに、5回表が出る確率を求めてみましょう。
表が出る回数が5回ですのでxは5, 10回投げるのでn=10, 正しいコインなのでp=0.5, 確率関数なので最後は0です。
前回、合計を求めるときに=SUM(A2:A31)と入力したように、どこかのセルに=BINOMDIST(5, 10, 0.5, 0)と入力します。
次に6回以上表が出る確率を求めてください。今回は講義中に正解を書けないので数字だけ書いておきます。0.376953…となれば正解、入力した式は正しいです。
しかし、BINOMDISTでは10,000回投げたときに5,000回以上表が出る確率を求めようとしてもnが大きすぎて出来ませんので、こんなときは正規分布で近似します。
正規分布N(μ,σ2)の密度関数f(x)はエクセルの関数NORMDIST(x,μ,σ,0)で求めることが出来ます。
一方分布関数P(X≦x)はNORMDIST(x,μ,σ,1)で求めることが出来ます。
練習:Xの確率分布が標準正規分布の時、P(X≦1)を求めなさい。0.841345…となれば正解です。
次に、二項分布B(n,p)でnが大きいときの正規分布による近似が正しいかどうか確かめるために、BINOMDISTでも計算できるn=100で確かめてみましょう。
正しいコインを100回投げたときに表が出た回数をXとします。
P(X≦50)とP(X≦55)をBINOMDISTを用いて求めてみましょう。それぞれ0.539795…と0.864373…となれば正解です。
Xの分布を正規分布で近似するために期待値μと分散σ2の値を求めてください。公式は教科書にまとめてあります。
P(X≦50)とP(X≦55)を正規分布による近似で求めてください。それぞれ0.5と0.841345…となれば正解です。
少し精度が低いですね。講義中に板書した、二項分布の棒グラフと正規分布の密度関数の曲線で囲まれた部分の面積を思い出してください。棒グラフには幅があるので、密度関数を積分するときにもその幅の分だけ広げた方が精度があがります。どのようにすれば精度が上がるか考えてください。0.539828…と0.864334…になります。
さて、n=10の場合はこのように確率関数、分布関数を計算することが出来ますし、正規近似では近似精度が悪いです。でもnが大きくなると確率関数を計算することが出来なくなる一方、正規近似の精度は良くなります。それで正規分布N(μ,σ^2)の分布関数で近似します。
とりあえずn=10の場合に対してどれくらい近似できるか見てみましょう。
まずx=0に対して
=NORMDIST(a6,d$18,sqrt(e$18),1)と入力します。3番目の引数には分散ではなく
その平方根である標準偏差を代入します。以下の画像では=NORMDIST(a6,d$18,sqrt(e$18),true)と書いています。本当はtrue, falseなのですが、面倒なのでtrueの代わりに1、falseの代わりに0を書いても同じ計算結果が得られるようになっています。
それをF16セルまでコピーします。
さて近似してみましたが、C列とF列を見比べてみると値はかなり異なっています。
今回の講義の最初の方に書いた、棒グラフの幅が原因です。特にn=10だと棒も11個しかないので、その幅の影響が大きいです。
まずx=0に対して
=NORMDIST(a6+0.5,d$18,sqrt(e$18),1)と入力します。つまり幅の分0.5増やしています。
それをG16セルまでコピーします。
如何ですか。小数点以下2桁までは合いました。まだ近似としては不十分ですが、
そもそもn=10ならば直接確率関数を計算した方が話は早いです。
ここまで出来た人は、このSheet1をSheet2へコピーして、n=100の場合について
同じ計算をして見ましょう。これ以上あまりnを大きくすると、確率関数が計算出来なくなったり、出来たとしても確率関数を足すときの誤差が大きくなってしまいます。
Sub 練習()
End Sub
となり、カーソルは真ん中の空行にあります。そこに
cells(1,1)=1
と書いてください。これは上から1番目、左から1番目のセルに1を代入しなさい、という命令です。
プログラムを書いたら、上の右向きの三角をクリックしてください。するとどのプログラムを実行するか聞かれます。今は一つしか書いていないのでそのまま実行を押すと、A1セルに1が代入されます。Excelの画面で「マクロ」を押しても実行できます。
このプログラムを保存する時の注意。
左上の「上書き保存」を押すと、どの形式でファイルを保存するか尋ねられます。ファイルの種類として「Excelマクロ有効ブック」を指定してください。
さてこれでプログラムを書いて実行できるようになったので、具体的なプログラミングを始めましょう。
先ほどのCells(1,1)=1とEnd Subの間に
cells(2,1)=2
cells(3,1)=3
と書き込んで実行してみてください。
あるいは
cells(1,2)=-2
cells(1,3)=-3
と書いて、cellsの二つの添え字がExcelのどのセルに対応しているか確認しましょう。
10個続けて書くには、
cells(1,1)=1
cells(2,1)=2
cells(3,1)=3
cells(4,1)=4
cells(5,1)=5
cells(6,1)=6
cells(7,1)=7
cells(8,1)=8
cells(9,1)=9
cells(10,1)=10
と書くのは大変なので、これをFor文を使って
for i=1 to 10
cells(i,1)=i
next i
と書きます。
次に下のような九九の表ができるようにプログラムを考えて見ましょう。
推定方法の良し悪しを比較するために、シミュレーションしてみましょう。
Sub tenchi()
Dim A As Range, B As Range
Set A = Range(Cells(1, 1), Cells(2, 2))
Set B = Range(Cells(4, 1), Cells(5, 2))
Cells(1, 1) = 11
Cells(1, 2) = 12
Cells(2, 1) = 21
Cells(2, 2) = 22
B = Application.Transpose(A)
End Sub
これでCells(1, 1)からCells(2, 2)までに入力した行列の転置がCells(4, 1)からCells(5, 2)までに出力されます。
Sub matrix()
Dim A As Range, B As Range, C As Range, D As Range
Cells(1, 1) = 11
Cells(1, 2) = 12
Cells(2, 1) = 21
Cells(2, 2) = 22
Set A = Range(Cells(1, 1), Cells(2, 2))
Set B = Range(Cells(4, 1), Cells(5, 2))
Set C = Range(Cells(7, 1), Cells(8, 2))
Set D = Range(Cells(10, 1), Cells(11, 2))
B = Application.Transpose(A)
C = Application.MInverse(A)
D = Application.MMult(A, C)
Cells(13, 1) = Application.MDeterm(D)
End Sub
となります。誤差が出るので0の筈なのに0にならない成分があるかもしれません。
正解はx=11, y=-6, z=5です。
全然関係ないですがカネゴンの78ちゃんねるまとめブログ、証拠を押さえられる | デジタルマガジン面白かったです。
そこで掛け算を足し算に変換するために、対数を考えることにします。すると
今日の株価 今日の株価 昨日の株価
－－－－－－＝－－－－－×－－－－－－
一昨日の株価 昨日の株価 一昨日の株価
は
Log(今日の株価)－Log(一昨日の株価)＝（Log(今日の株価)－Log(昨日の株価)）＋（Log(昨日の株価)－Log(一昨日の株価)）
という、足し算、引き算の式になります。
これを先ほど計算したF列の隣、G列に書きます。
対数の底は10でもeでも上記の計算は成り立ちますが、底がe、つまり自然対数の方が
log(x)≒x-1
という近似式が成り立つので、
log(今日の株価÷昨日の株価)≒今日の株価÷昨日の株価－１≒(今日の株価－昨日の株価)÷昨日の株価
なので先ほど求めた変化率と似た値になります。
計算
G2セルに
=LN(E2)-LN(E3)
と書いても計算できますが、先ほど同様VBAを使って計算して下さい。但し自然対数はExcelのセルに書くときはLnですが、VBAで書くときにはLogです。
Windows7のスクリーンショット
バージョンが6.1と書かれています。
一つ疑問があるのですが、マイクロソフトユーザー向け特別優待プログラムって、そもそもWindows用のこのソフトをマイクロソフトユーザー以外が買うことあるのでしょうか？
