-
Notifications
You must be signed in to change notification settings - Fork 5
/
16.html
811 lines (728 loc) · 31.2 KB
/
16.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
<!doctype html>
<html lang="ja">
<head>
<meta charset="utf-8">
<title>プログラマの為の数学勉強会</title>
<!-- For reveal.js -->
<link rel="stylesheet" href="lib/reveal/css/reveal.css">
<link rel="stylesheet" href="lib/reveal/css/theme/night.css">
<link rel="stylesheet" href="lib/reveal/lib/css/ir_black.css">
<!-- For Graphics -->
<link rel="stylesheet" href="css/graphics.css">
<style>
.reveal .chapter-title {
margin-top: 3em;
}
.reveal {
font-size: 32px;
line-height: 1.4em;
}
.reveal .slides {
text-align: left;
}
.reveal section img {
border: none;
background: 0;
margin-left: 1em;
margin-right: 1em;
box-shadow: none;
}
.reveal strong {
color: yellow;
}
.reveal sup {
font-size: 40%;
}
.reveal table {
margin-top: 0.5em;
margin-bottom: 0.5em;
border: 2px solid lightblue;
}
.reveal pre {
font-size: 0.7em;
}
.reveal pre code {
max-height: 600px;
}
.reveal .note {
font-size: 50%;
}
.reveal .controls div.navigate-up,
.reveal .controls div.navigate-down {
display: none;
}
.reveal .block {
border: solid 2px;
position: relative;
border-radius: 8px;
margin-top: 0.8em;
margin-bottom: 0.8em;
padding: 1em 0.8em 1em 0.8em;
}
.reveal .block:after {
content: "";
display: block;
clear: both;
height: 1px;
overflow: hidden;
}
.reveal .answer {
color: #111111;
}
.reveal .block h4 {
position: absolute;
top: -0.5em;
margin: 0 auto;
background: #111111;
font-weight: bold;
}
</style>
<!-- Setup libraries for RequireJS-->
<script src="lib/require.js"></script>
<script>
requirejs.config({
baseUrl: "js",
paths: {
d3: "../lib/d3/d3.v3.min",
numeric: "../lib/numeric-1.2.6",
MathJax: "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML"
},
shim: {
d3: { exports: "d3" },
numeric: { exports: "numeric" },
MathJax: { exports: "MathJax" }
}
});
</script>
<!-- Initialize MathJax -->
<script type="text/x-mathjax-config">
require(["MathJax"], function (MathJax){
MathJax.Hub.Register.StartupHook("AsciiMath Jax Config", function () {
var AM = MathJax.InputJax.AsciiMath.AM;
AM.symbols.push(
{input:"mathbi",tag:"mstyle",atname:"mathvariant",atval:"bold-italic",
output:"mathbi",tex:null,ttype:AM.TOKEN.UNARY}
);
});
MathJax.Hub.Config({
showProcessingMessages: false,
skipStartupTypeset: false,
tex2jax: {
inlineMath: [ ["\\(","\\)"] ],
displayMath: [ ["\\[","\\]"] ]
}
});
});
</script>
<script>
</script>
</head>
<body>
<div class="reveal">
<div class="slides">
<section style="text-align: center">
<h1> プログラマの為の<br>数学勉強会<br>第16回</h1>
<span>
(於)ワークスアプリケーションズ<br>
中村晃一<br>
2014年1月16日
</span>
</section>
<section>
<h2>謝辞</h2>
<p>
この会の企画・会場設備の提供をして頂きました<br>
㈱ ワークスアプリケーションズ様<br>
にこの場をお借りして御礼申し上げます。
</p>
</section>
<section>
<h2> この資料について </h2>
<p>
<ul>
<li> <a href="http://nineties.github.com/math-seminar">
http://nineties.github.com/math-seminar
</a>に置いてあります。 </li>
<li> SVGに対応したブラウザで見て下さい。主要なブラウザで古いバージョンでなければ大丈夫だと思います。</li>
<li> 内容の誤り、プログラムのバグは<a href="http://twitter.com/9_ties">@9_ties</a>かkoichi.nakamur AT gmail.comまでご連絡下さい。</li>
<li> サンプルプログラムはPython及びMaximaで記述しています。 </li>
</ul>
</p>
</section>
<section>
<h2 class="chapter-title"> 適合度検定・独立性検定 </h2>
</section>
<section>
<p>
前回は, 母集団分布の形が既知の場合に母平均や母分散の値, それらに差があるかといった事を検定する方法を紹介しました. これらは母数というパラメータを考えるので <strong>パラメトリック</strong> 手法と呼ばれます.
</p>
<p>
今回紹介するのは <strong> ノンパラメトリック </strong>手法と呼ばれる, 母集団分布の形を仮定しない検定法のうち度数分布に適用出来る手法を紹介します.
</p>
</section>
<section>
<p>
検定の基本的な流れは同じです. 復習すると以下のように行うのでした.
</p>
<ol>
<li> 帰無仮説\(H_0\)と対立仮説\(H_1\)を立てる. </li>
<li> 有意水準 \(\alpha\) を決める. </li>
<li> 検定に用いる統計量を定めて, \(H_0,H_1\) に基き棄却域を計算する. </li>
<li> 実現値(実際に得られたデータ)に基いて, \(H_0\) を棄却するか否かを決める. </li>
</ol>
</section>
<section>
<h2> \(\chi^2\) 適合度検定 </h2>
<p>
<strong> 度数分布 </strong>, つまり～が \(x\) 回生じたという形の確率分布に対しては \(\chi^2\) 適合度検定を利用する事が出来ます.
</p>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 定理 </h4>
<p>
ある度数分布に従う確率変数のカテゴリ毎の観測値が \(x_1,x_2,\ldots,x_m\) 回, 理論値が \(e_1,e_2,\ldots,e_m\) である場合
</p>
\[ \chi^2 = \frac{(x_1 - e_1)^2}{e_1} + \frac{(x_2-e_2)^2}{e_2} + \cdots + \frac{(x_m-e_m)^2}{e_m} \]
は近似的に自由度 \(m-1\) の \(\chi^2\) 分布に従う.
</div>
<p>
\(\chi^2\) の値が大きいほど理論値から外れているという事になりますので, この定理に基いて\(\chi^2\)分布を利用した片側検定を行えば良いです.
</p>
</section>
<section>
<div class="block" style="border-color:lightgreen;font-size:90%">
<h4 style="color:lightgreen"> 例題 </h4>
<p>
サイコロを600回投げた所, それぞれの出目の回数が以下のようになりました. このサイコロは一様分布に従っていると言えるでしょうか？
\[ \begin{array}{|c|c|c|c|c|c|} \hline
1 & 2 & 3 & 4 & 5 & 6 \\ \hline
112\text{回} & 85\text{回} & 120\text{回} & 101\text{回} & 89\text{回} & 93\text{回} \\ \hline
\end{array} \]
</p>
</div>
<p class="fragment" style="font-size:80%">
帰無仮説・対立仮説を
\[ \begin{aligned}
H_0&: \text{分布は一様である.} \\
H_1&: \text{分布は一様ではない.}
\end{aligned} \]
と仮定します.
</p>
</section>
<section style="font-size:80%">
<p> 仮説 \(H_0\) の元では, 各目は確率 \(1/6\) なので
\[ \begin{aligned}
\chi^2 &= \frac{(112 - 100)^2}{100} + \frac{(85-100)^2}{100} + \frac{(120-100)^2}{100} + \frac{(101-100)^2}{100} \\
& + \frac{(89-100)^2}{100} + \frac{(93-100)^2}{100} = 9.4
\end{aligned} \]
となります.
</p>
<p>
ここで, 自由度 \(5\) の \(\chi^2\) 分布の上側 5% 点を求めると
</p>
<pre><code class="python" style="max-height:400px">>>> from scipy import stats
>>> stats.chi2.ppf(0.95, 5)
11.070497693516351
</code></pre>
<p>
となるので, \(H_0\) を棄却する事は出来ません.
</p>
</section>
<section>
<div class="block" style="border-color:lightgreen;font-size:90%">
<h4 style="color:lightgreen"> 練習問題 </h4>
あるWebサービスへの1秒あたりのアクセス回数の測定を1000回行いました. これはポアソン分布 \(Po(5)\) に従っていると言えるでしょうか？
\[ \tiny{ \begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline
\text{0回} & \text{1回} & \text{2回} & \text{3回} & \text{4回} & \text{5回} & \text{6回} & \text{7回} & \text{8回} & \text{9回} & \text{10回} & \text{11回} \\ \hline
5 & 31 & 72 & 135 & 183 & 173 & 162 & 89 & 79 & 41 & 18 & 12 \\ \hline
\end{array}} \]
</div>
<p class="fragment" style="font-size:70%">
帰無仮説 \(H_0\) を母集団分布が \(Po(5)\) である事とします. 1秒間のアクセス回数が \(k\) 回である確率は \(5^k e^{-5}/k!\) であるので理論値は
\[ \tiny{ \begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline
\text{0回} & \text{1回} & \text{2回} & \text{3回} & \text{4回} & \text{5回} & \text{6回} & \text{7回} & \text{8回} & \text{9回} & \text{10回} & \text{11回} & \text{12回} & \text{13回} & \text{14回} & \cdots \\ \hline
6.7 & 33.7 & 84.2 & 140.4 & 175.5 & 175.5 & 146.2 & 104.4 & 65.3 & 36.3 & 18.1 & 8.2 & 3.4 & 1.3 & 0 & \cdots\\ \hline
\end{array}} \]
となります. これから \(\chi^2 = 17.0\) となります. 自由度 \(13\) の\(\chi^2\) 分布の上側 5% 点は \(22.4\) なので, \(H_0\) を棄却する事は出来ません.
</p>
</section>
<section>
<h2> 独立性検定 </h2>
<p>
適合度検定を応用して <strong>独立性検定 </strong> を行う事が出来ます.
</p>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 定理 </h4>
<p>
2つの独立な度数分布に従う確率変数 \(X(\text{カテゴリ数 $m$})\),\(Y(\text{カテゴリ数 $n$})\) について
\[ \chi^2 = \sum \frac{(\text{観測値}-\text{理論値})^2}{\text{理論値}} \]
</p>
は近似的に自由度 \((m-1)\times(n-1)\) の \(\chi^2\) 分布に従う.
</div>
</section>
<section>
<div class="block" style="border-color:lightgreen;font-size:90%">
<h4 style="color:lightgreen"> 例題 </h4>
<p>
新薬の効果を見るために新薬と古い薬を適当な人数ずつ処方したところ, 以下の結果となった. 新薬に効果はあったと言えるか？
</p>
\[ \begin{array}{|c|cc|c|} \hline
& 効果あり & 効果なし & \text{計} \\ \hline
\text{新しい薬} & 79 & 41 & 120 \\
\text{古い薬} & 43 & 37 & 80 \\ \hline
\text{計} & 122 & 78 & 200 \\ \hline
\end{array} \]
</div>
</section>
<section style="font-size:70%">
<p>
新薬でも古い薬でも効果が変わらないならば, \((\text{効果あり}):(\text{効果なし})\) の比率は一定です. そこで, \((\text{効果あり}):(\text{効果なし}) = 122:78 = 61:39\) であると仮定するならば, 理論値は
\[ \begin{array}{|c|cc|c|} \hline
& 効果あり & 効果なし & \text{計} \\ \hline
\text{新しい薬} & 73.2 & 46.8 & 120 \\
\text{古い薬} & 48.8 & 31.2 & 80 \\ \hline
\text{計} & 122 & 78 & 200 \\ \hline
\end{array} \]
となるはずです. よって
\[ \chi^2 = \frac{(79-73.2)^2}{73.2}+\frac{(41-46.8)^2}{46.8}+\frac{(43-48.8)^2}{48.8}+\frac{(37-31.2)^2}{31.2} = 2.9 \]
となります. よって自由度 \((2-1)\times(2-1)=1\) の \(\chi^2\) 分布の上側 5% 点は \(3.8\) であるので, 新薬に効果があったとは言えません.
</p>
</section>
<section>
<h2 class="chapter-title"> 回帰分析 </h2>
</section>
<section>
<p>
ある結果が幾つかの要因からどのように説明されるかを推定する手法を <strong>回帰分析</strong> といいます. 例えば以下の様な問題を考えます.
</p>
<div class="block fragment" style="border-color:pink;font-size:90%">
あるWebサービスでのユーザーの平均滞在時間は年齢・性別・職業などの要因によってどのように決定されるだろうか？
</div>
<p class="fragment">
数学的には, 確率変数 \(Y\) と確率変数 \(X_1,X_2,\ldots,X_n\) について
\[ Y=f(X_1,X_2,\ldots,X_n) \]
を満たすような \(f\) を推定する事となります. この \(f\) を \(Y\) の <strong> モデル </strong> と呼びます.
</p>
</section>
<section>
<h2> 最小二乗法(1変数) </h2>
<p>
観測値 \(Y\) と理論値 \(f(X)\) の誤差 \(Y-f(X)\) が正規分布に従うと仮定して最尤推定で \(f(X)\) を決定する手法を <strong> 最小二乗法 </strong> と呼びます.
</p>
</section>
<section style="font-size:80%">
<p>
誤差 \(E=Y-f(X)\) が正規分布 \( N(0,\sigma^2) \) に従うと仮定します. すると, \(X=x_i\) と値が確定した場合 \(Y\) は \(N(f(x_i),\sigma^2)\) に従います. よって, データ列 \((x_1,y_1),(x_2,y_2),\ldots,(x_n,y_n)\) が得られたならば, 尤度関数は
\[ \begin{aligned}
L(f) &= \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(y_i-f(x_i))^2}{2\sigma^2}\right) \\
&\propto \exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-f(x_i))^2\right)
\end{aligned} \]
となります. よって \(L(f)\) が最大になるのは <strong> 残差平方和 </strong>
\[ \sum_{i=1}^n(y_i-f(x_i))^2 \]
が最小となる場合となります.
</p>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 最小二乗法(1変数・母分散が一定) </h4>
<p>
誤差の母分散が一定の場合, データ列 \((x_1,y_1),(x_2,y_2),\ldots,(x_n,y_n)\) に対して残差平方和
\[ \sum_{i=1}^n\left\{ y_i-f(x_i) \right\}^2 \]
を最小とする \(f\) としてモデルを決定する.
</p>
</div>
</section>
<section>
<div class="block" style="border-color:pink;font-size:80%">
<h4 style="color:pink"> 最小二乗法:一次関数の場合 </h4>
<p>
モデルが \(y=ax+b\) という形の場合, データ列 \((x_1,y_1),(x_2,y_2),\ldots,(x_n,y_n)\) に対して残差平方和を最小とする\(a,b\) は
\[ a = \frac{n\sum x_iy_i - \sum x_i \sum y_i}{n\sum x_i^2 - \left( \sum x_i\right)^2}, b = \frac{\sum x_i^2 \sum y_i - \sum x_iy_i \sum x_i}{n\sum x_i^2 - \left( \sum x_i \right)^2} \]
である.
</p>
</div>
</section>
<section style="font-size:70%">
<p>
【証明】<br>
\[ E=\sum_{i=1}^n(y_i-ax_i-b)^2 \]
を \(a,b\) で偏微分すれば
\[ \begin{aligned}
\frac{\partial E}{\partial a} &= -\sum 2x_i(y_i-ax_i-b)=2(\sum x_i^2)a+2(\sum x_i) b - 2\sum x_iy_i \\
\frac{\partial E}{\partial b} &= -\sum 2(y_i-ax_i-b) = 2(\sum x_i)a + 2nb - 2\sum y_i
\end{aligned} \]
となるので, \(E\) が最小となるのは
\[ \begin{aligned}
&\frac{\partial E}{\partial a}=\frac{\partial E}{\partial b} = 0 \\
\Leftrightarrow &
a = \frac{n\sum x_iy_i - \sum x_i \sum y_i}{n\sum x_i^2 - \left( \sum x_i\right)^2}, b = \frac{\sum x_i^2 \sum y_i - \sum x_iy_i \sum x_i}{n\sum x_i^2 - \left( \sum x_i \right)^2}
\end{aligned} \]
の時である. <span style="float:right">□</span>
</p>
</section>
<section>
<div class="block" style="border-color:lightgreen;font-size:80%">
<h4 style="color:lightgreen"> 練習問題 </h4>
<p>
以下のデータ列に対して \(y=ax+b\) というモデルを仮定して, 最小二乗法によって \(a,b\) を決定して下さい.
</p>
\[ \small{\begin{array}{|c|cccccccccc|} \hline
x & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\ \hline
y & 1.6 & 2.9 & 7.8 & 11.2 & 11.7 & 14.2 & 16.0 & 20.0 & 17.7 & 18.9 \\ \hline
\end{array}} \]
<div align="center">
<img width="400px" src="fig/ex16-1.png">
</div>
</div>
</section>
<section>
<p>
(答え) \( a=2.0, b=3.0 \)
</p>
<div align="center">
<img width="600px" src="fig/ex16-1-2.png">
</div>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 最小二乗法(1変数・母分散が既知) </h4>
<p>
データ列 \((x_1,y_1),(x_2,y_2),\ldots,(x_n,y_n)\) に対して, \(i\) 番目のデータの誤差の母分散が \(\sigma_i^2\) であるならば, 残差平方和
\[ \sum_{i=1}^n\frac{\left\{ y_i-f(x_i) \right\}^2}{\sigma_i^2} \]
を最小とする \(f\) としてモデルを決定する.
</p>
</div>
<p> 【練習問題】\(y=ax+b\) というモデルの場合に \(a,b\) の値を得る公式を求めて下さい. </p>
</section>
<section>
<p>
数値計算の手法を利用すれば, 一次式でない場合にも最小二乗法を適用する事が可能です.
</p>
</section>
<section>
<div class="block" style="border-color:lightgreen;font-size:80%">
<h4 style="color:lightgreen"> 練習問題 </h4>
<p>
以下のデータ列に対して \(y=\sin ax \) というモデルを仮定して, 最小二乗法によって \(a\) を決定して下さい.
</p>
\[ \small{\begin{array}{|c|cccccccccc|} \hline
x & 0.0 & 0.5 & 1.0 & 1.8 & 2.1 & 2.6 & 3.1 & 3.7 & 4.2 & 4.7 \\ \hline
y & 0.0 & 0.2 & 0.5 & 0.7 & 0.9 & 1.1 & 1.1 & 1.1 & 0.9 & 0.6 \\ \hline
\end{array}} \]
<div align="center">
<img width="400px" src="fig/ex16-2.png">
</div>
</div>
</section>
<section style="font-size:80%">
<p> 第4回にやったニュートン法を復習しておきます. </p>
<div class="block" style="border-color:pink;font-size:80%">
<h4 style="color:pink"> ニュートン法 </h4>
<p>
方程式\(f(x)=0\)の解を\(\alpha\)とする。
</p>
<p>
\(x=\alpha\)を含む適当な区間で\(f'(x)\neq 0\),\(f(x),f'(x),f''(x)\)が連続ならば,この区間内に初期値\(x_0\)を取れば数列
\[ x_{n+1} = x_{n} - \frac{f(x_n)}{f'(x_{n})} \]
が\(\alpha\)に収束する。
</p>
</div>
<div align="center"> <img width="300px" src="fig/newton-method.png"> </div>
</section>
<section style="font-size:80%">
<p>
\[ E = \sum (y_i-\sin ax_i)^2 \]
を \(a\) について微分すれば
\[ \frac{\mathrm{d} E}{\mathrm{d} a} = -2\sum x_i\cos ax_i (y_i-\sin ax_i) = -2 \sum x_iy_i\cos ax_i+\sum x_i\sin 2ax_i \]
となるので, 方程式
\[ \sum x_i \sin 2ax_i - 2\sum x_iy_i\cos ax_i = 0\]
を解けば良いです. この左辺を \(f(a)\) とおくと
\[ f'(a) = \sum 2x_i^2\cos 2ax_i + 2\sum x_i^2y_i\sin ax_i \]
となるので, ニュートン法
\[ a_{k+1} = a_k - \frac{f(a_k)}{f'(a_k)}=a_k - \frac{\sum x_i \sin 2a_kx_i - 2\sum x_iy_i\cos a_kx_i}{\sum 2x_i^2\cos 2a_kx_i + 2\sum x_i^2y_i\sin a_kx_i} \]
を利用する事が出来ます. 初期値はグラフを見て \(a_0 = 0.4\)くらいにすれば良いでしょう.
</p>
</section>
<section>
<p>
以下がコーディング例です.
</p>
<pre><code class="python" style="font-size:90%;max-height:400px">import numpy as np
x = np.array([0.0, 0.5, 1.0, 1.8, 2.1, 2.6, 3.1, 3.7, 4.2, 4.7])
y = np.array([0.0, 0.2, 0.5, 0.7, 0.9, 1.1, 1.1, 1.1, 0.9, 0.6])
a = 0.4
while True:
next_a = a - (np.sum( x*np.sin(2*a*x) ) - 2 * np.sum( x*y*np.cos(a*x) )) / (np.sum( 2*x*x*np.cos(2*a*x) ) + 2 * np.sum( x*x*y*np.sin(a*x) ))
if abs((next_a-a)/a) < 1.0e-6:
break
a = next_a
print a
</code></pre>
</section>
<section>
<p>
計算結果は \(a=0.5\) となります.
</p>
<div align="center">
<img width="600px" src="fig/ex16-2-2.png">
</div>
</section>
<section>
<div class="block" style="border-color:lightgreen;font-size:80%">
<h4 style="color:lightgreen"> 練習問題 </h4>
<p>
以下のデータ列に対して \(y=\sin ax \) というモデルを仮定して, 最小二乗法によって \(a\) を決定して下さい.
</p>
\[ \small{\begin{array}{|c|cccccccccc|} \hline
x & 0.0 & 0.5 & 1.0 & 1.8 & 2.1 & 2.6 & 3.1 & 3.7 & 4.2 & 4.7 \\ \hline
y & 0.0 & 0.2 & 0.5 & 0.7 & 0.9 & 1.1 & 1.1 & 1.1 & 0.9 & 0.6 \\ \hline
\end{array}} \]
<div align="center">
<img width="400px" src="fig/ex16-2.png">
</div>
</div>
</section>
<section>
<div class="block" style="border-color:pink;font-size:90%">
<h4 style="color:pink"> 最小二乗法(線型回帰・母分散が一定) </h4>
<p>
モデル \(y=f(x)\) が既知の関数 \(g_1(x),g_2(x),\ldots,g_m(x)\)とパラメータ \(a_1,a_2,\ldots,a_m\) によって
\[ f(x)=a_1g_1(x) + a_2g_2(x) + \cdots + a_mg_m(x) \]
と表される場合,
\[\begin{aligned}
G&=(g_j(x_i))_{1\leq i \leq n,1\leq j\leq m}\\
\mathbf{a}&=(a_1,a_2,\ldots,a_m)^T\\
\mathbf{y}&=(y_1,y_2,\ldots,y_n)^T
\end{aligned} \]
とおけば <strong> 正規方程式 </strong>
\[ \color{yellow}{ G^TG\mathbf{a}=G^T\mathbf{y} } \]
の解が最小二乗法による \(a_1,a_2,\ldots,a_m\) の推定値を与える.
</p>
</div>
</section>
<section>
<p>
例えば,
\[ f(x)=a_nx^n + a_{n-1}x^{n-1} + \cdots + a_1x + a_0 \]
であるとか,
\[ f(x) = a_n\sin nx + a_{n-1}\sin (n-1)x + \cdots + a_1\sin x + a_0 \]
などの形のモデルであれば今の方法を使う事が出来ます.
\[ f(x)= e^{ax} \]
の様に, 関数の中にパラメータが入っている場合には使えないので数値解法を利用したり, モデルの法を近似するなどの工夫を行って下さい.
</p>
</section>
<section>
<p>
証明の前に<a href="http://nineties.github.com/math-seminar/8.html#/62">第8回</a>にやった<strong>ベクトルでの微分法</strong>を復習しておきます.
</p>
<div class="block" style="border-color:pink;font-size:80%">
<h4 style="color:pink"> ベクトルでの微分法 </h4>
<p>
多変数スカラー値関数\(f,g\),定数\(k,l\)に関して
\[ \frac{\mathrm{d}}{\mathrm{d}\mathbf{x}}(kf+lg) = k\frac{\mathrm{d}f}{\mathrm{d}\mathbf{x}} + l\frac{\mathrm{d}g}{\mathrm{d}\mathbf{x}} \]
また\(c\)を定数,\(\mathbf{a}\)を定ベクトル,\(A\)を正方行列とすると
\[ \begin{aligned}
&\frac{\mathrm{d}}{\mathrm{d}\mathbf{x}}c = \mathbf{0} \\
&\frac{\mathrm{d}}{\mathrm{d}\mathbf{x}}\mathbf{a}^T\mathbf{x} = \frac{\mathrm{d}}{\mathrm{d}\mathbf{x}}\mathbf{x}^T\mathbf{a} = \mathbf{a} \\
&\frac{\mathrm{d}}{\mathrm{d}\mathbf{x}}||\mathbf{x}||^2 = 2\mathbf{x} \\
&\frac{\mathrm{d}}{\mathrm{d}\mathbf{x}}\mathbf{x}^TA\mathbf{x} = (A+A^T)\mathbf{x}
\end{aligned} \]
が成り立つ。
</p>
</div>
</section>
<section>
<p>
【証明】 <br>
残差平方和
\[ E=\sum(y_i-a_1g_1(x_i)-a_2g_2(x_i)-\cdots-a_mg_m(x_i))^2 \]
はベクトル \(\mathbf{y}-G\mathbf{a}\) のノルムの二乗であるから
\[ \begin{aligned}
E&=(G\mathbf{a}-\mathbf{y})^T(G\mathbf{a}-\mathbf{y}) = (\mathbf{a}^TG^T-\mathbf{y}^T)(G\mathbf{a}-\mathbf{y}) \\
&= \mathbf{a}^TG^TG\mathbf{a}-\mathbf{a}^TG^T\mathbf{y}-\mathbf{y}^TG\mathbf{a}+\mathbf{y}^T\mathbf{y}
\end{aligned} \]
となるので,
\[ \frac{\mathrm{d} E}{\mathrm{d} \mathbf{a}} = (G^TG+(G^TG)^T)\mathbf{a}-2G^T\mathbf{y} = 2(G^TG\mathbf{a}-G^T\mathbf{y}) \]
である. よって \(E\) が最小値を取る点で
\[ G^TG\mathbf{a} = G^T\mathbf{y} \]
が成立する. <span style="float:right">□</span>
</p>
</section>
<section>
<div class="block" style="border-color:lightgreen;font-size:80%">
<h4 style="color:lightgreen"> 練習問題 </h4>
<p>
以下のデータに対して \(y=ax^2+bx+c\) というモデルを最小二乗法によって当てはめて下さい.
</p>
<pre><code class="python" style="font-size:90%;max-height:400px">#x y
0.0 3.7
0.5 1.2
1.0 -3.2
1.5 -1.0
2.0 -4.1
2.5 -3.3
3.0 -3.5
3.5 -2.4
4.0 -1.8
4.5 1.7
5.0 2.7
5.5 5.5
6.0 8.5
6.5 11.4
7.0 17.1
7.5 22.4
8.0 26.4
8.5 33.2
9.0 39.2
9.5 46.2
</code></pre>
</div>
</section>
<section>
<p>
モデルが \(y=ax^2+bx+c\) という形なので,
\[ g_1(x) = x^2, g_2(x) = x, g_3(x) = 1 \]
とおきます. すると
\[ G = \begin{pmatrix}
x_1^2 & x_1 & 1 \\
x_2^2 & x_2 & 1 \\
\vdots & \vdots & \vdots \\
x_n^2 & x_n & 1
\end{pmatrix},
\mathbf{a} = \begin{pmatrix}
a \\ b \\ c
\end{pmatrix},
\mathbf{y} = \begin{pmatrix}
y_1 \\ y_2 \\ \vdots \\ y_n
\end{pmatrix} \]
となります. これに基いて正規方程式を導出し, それを解きます.
</p>
</section>
<section>
<p>
以下がコーディング例です.
</p>
<pre><code class="python" style="font-size:90%;max-height:400px">import numpy as np
from scipy import linalg as LA
N = 20
x = np.array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5])
y = np.array([3.7, 1.2, -3.2, -1.0, -4.1, -3.3, -3.5, -2.4, -1.8, 1.7, 2.7, 5.5, 8.5, 11.4, 17.1, 22.4, 26.4, 33.2, 39.2, 46.2])
G = np.array([x*x, x, np.ones(N)]).T
print LA.solve(G.T.dot(G), G.T.dot(y))
</code></pre>
</section>
<section>
<p> 結果は
\[ a= 1.0, b= -5.2, c=3.1 \]
となります.
</p>
<div align="center">
<img width="600px" src="fig/ex16-3-2.png">
</div>
</section>
<section>
<p>
測定誤差の母分散が一定ではないが既知という場合にも, 一変数の場合と同様に解くことが出来ます. 方程式を導出してみてください.
</p>
</section>
<section>
<h2> 今回はここで終わります。 </h2>
<p>
次回も回帰分析の続きを行います. モデルの当てはまりの良さを判定や, 多変量回帰分析の計算法, 過学習といった問題の紹介を行います.
</p>
</section>
</div>
</div>
<script src="lib/reveal/lib/js/head.min.js"></script>
<script src="lib/reveal/js/reveal.js"></script>
<script>
Reveal.initialize({
width: 960,
height: 640,
controls: true,
progress: false,
history: true,
overview: false,
touch: true,
center: false,
rollingLinks: false,
transition: "page",
transitionSpeed: "default",
// When scale != 1, positions of mouse events will be incorrect.
minScale: 1.0,
maxScale: 1.0,
dependencies: [
{ src: "lib/reveal/lib/js/classList.js", condition: function() { return !document.body.classList; } },
{ src: "lib/reveal/plugin/highlight/highlight.js", async: true, callback: function() { hljs.initHighlightingOnLoad(); } },
{ src: "lib/reveal/plugin/zoom-js/zoom.js", async: true, condition: function() { return !!document.body.classList; } },
{ src: "lib/reveal/plugin/notes/notes.js", async: true, condition: function() { return !!document.body.classList; } }
]
});
// register event listeners
require(["MathJax"], function (MathJax) {
// Delay typesetting of slides
function typeset (idx) {
for (var i = idx - 2, n = idx + 2; i <= n; i++) {
var slide = Reveal.getSlide(i);
if (!slide) continue;
if (!slide.typeset) {
MathJax.Hub.Typeset(slide);
slide.typeset = true;
}
}
}
function initializeGraphics (idx) {
for (var i = idx - 2, n = idx + 2; i <= n; i++) {
var slide = Reveal.getSlide(i);
if (!slide) continue;
var graphics = slide.getAttribute("graphics");
if (graphics && !slide.initialized) {
slide.initialized = true;
(function () {
var p = slide;
require([graphics], function(g) {
if (g.initialize) g.initialize(p);
});
})();
}
}
}
function start (slide) {
var graphics = slide.getAttribute("graphics");
if (graphics) {
require([graphics], function(g) { if (g.start) g.start(slide); });
}
}
function stop (slide) {
var graphics = slide.getAttribute("graphics");
if (graphics) {
require([graphics], function(g) { if (g.stop) g.stop(slide); });
}
}
function simpleEvent (type) {
var event = document.createEvent("HTMLEvents");
event.initEvent(type, true, true);
return event;
}
Reveal.addEventListener("slidechanged", function (event) {
typeset(event.indexh);
initializeGraphics(event.indexh);
start(event.currentSlide);
stop(event.previousSlide);
});
Reveal.addEventListener("fragmentshown", function (event) {
var slide = Reveal.getCurrentSlide();
var graphics = slide.getAttribute("graphics");
if (graphics) {
require([graphics], function(g) { if (g.proceed) g.proceed(slide); });
}
});
console.log(Reveal.getIndices().h);
typeset(Reveal.getIndices().h);
initializeGraphics(Reveal.getIndices().h);
start(Reveal.getCurrentSlide());
});
</script>
</body>
</htl>