傾きの検定

温暖化しているかどうか、傾きの検定をしました。

二つの母集団の期待値に関する検定

2つの母分散の比に関する検定と推定

二つの母集団の母分散に関する検定と推定を説明しました。

季節成分

季節成分の計算方法について説明しました。

オートフォーカスの合いやすさ比較

デジカメで暗いところとか、飛行機雲のように細いもの、小さいものを撮ると、ピントが合わないことがありました。
デジカメのカタログを見ると、明るい場所でのオートフォーカスの速さを競っていますが、ピントを合わせやすい被写体に対して合わせる時間が例えば0.1秒から0.01秒まで10分の1になることよりも、実際に使うときは速さよりも暗い場所や小さい被写体にもちゃんとピントが合うことの方がもっと大事です。
デジカメのオートフォーカスは、ピントが合いやすい順に「位相差AF」「像面位相差AF」「コントラストAF」の三種類があります。
位相差AFが最もピントが合いやすいのですが、一眼レフだけしか使えません。一眼レフでも動画を撮るときには位相差AFを使えないのでコントラストAFになります。
コンデジやミラーレス一眼、それにスマホの殆どはコントラストAFで、一部の機種だけが像面位相差AFとコントラストAFの併用です。
一眼レフは大きくて旅行には邪魔になることがあるので、像面位相差AFやコントラストAFではどれくらいピントが合うかどうかを、ピントが合いにくいイルミネーションで試してみました。

まずは普通のコントラストAFのCanon PowerShot G7X

イルミネーションはコントラストAFではピントを合わせづらいです。

次は、コントラストAFの中では一番ピント合わせが速いと言われているPanasonic Lumix LX100

ピントを合わせやすい場面では速いのですが、イルミネーションには合わせにくいです。

像面位相差AFのSony α6000

ピントが正確ですし、合わせるのも速いです。
ただ、ブログに載せるために縮小すると、ピンボケの方がイルミネーションの光が大きく写って綺麗ですね。ピントがちゃんとあってるα6000の写真の方が寂しく感じます。
でも画像をクリックして拡大すると、ピンボケ写真はやっぱり残念です。
α6000でも、一眼レフと比べたらピントが合いにくいので、α6000を使うのは一眼レフが大きくて持ち運びにくい時だけです。
なのでα6000には一番小さなレンズを着けています。大きなレンズを着けるくらいなら最初から一眼レフを使うので。

最後に位相差AFの一眼レフD7100

持って行けるならこれが一番いいです。でも大きくて重たいです。

スマートフォンのカメラは、小さいので画質はコンデジに負けるはずですが、小さいのでピントが合う幅が広く、ピンボケになりにくいです。
どんなに画質の良いカメラでも、ピンボケになったら意味ないので、像面位相差AFの使えないコンデジ使うくらいなら最初からスマホで撮ったほうが良いと思います。
ズームで拡大すると、スマホはボケてきますので、コンデジの方が向いています。

母平均の推定と検定

1つの母集団の平均について、検定と推定を説明しました。

季節成分

移動平均との差に注目して、季節成分を計算しました。

母分散の推定と検定

今日から統計解析の実例です。今日は1標本母分散の推定と検定を説明しました。

時系列解析

前回までに説明した方法で取得したデータを用いて、時系列解析の方法を説明します。

SIMフリー機でドコモメールを使う

ドコモでデータ通信するにはspモード(324円/月)かmoperaU(540円/月)のどちらかを契約する必要があります。spモードの方が安いのですが、これまではspモードを使えるのはドコモの端末だけでした。
しかしドコモ：spモードを他社端末でも利用可能に – SIMロック解除義務化でということでHuawae Ascend Mate7など日本国内で販売されているSIMフリースマホでも使えるようになりました。
手続きはdocomo IDさえ作っていればMy docomoから簡単に出来ます。作っていなくてもSMSを受信できる端末(殆どのスマホで出来ます)を持っていれば作れます。
spモードを契約したらAPNを設定すればデータ通信出来ます。APNは spmode.ne.jp 認証タイプは CHAP です。IDやパスワードは不要です。

spモードを契約したのでデータ通信だけでなく @docomo.ne.jp のドコモメールも使えるようになりました。
最近は少なくなりましたが、キャリアのメールアドレスしか登録に使えないことがありましたので @docomo.ne.jp のメールアドレスが使えるなら使おうと思います。
ドコモのスマホを持っていれば簡単なのですが、SIMフリー機しか持っていないと、使うまでの準備が大変でした。
ドコモの説明を見ても、ドコモのスマホを持っている人向けにしか書かれていないので、SIMフリー機で使おうとして何故使えないのか、どの設定が足りないのか探すのに一苦労。
まるでアドベンチャーゲームでどのフラグが足りなくて先に進めないのかを考えるようなもので、東京からの帰りの新幹線の中で3時間あれこれ試してみてやっと使えるようになりました。
一番参考になったのはiモードメールをスマホで送受信する（iPhone偽装による設定方法）でしたが、UserAgent偽装は必要ありませんでした。

手順ですが、まずdocomo IDを作ってspモードを契約します。どちらもMy docomoで出来ます。
spモードを契約したら、spモードでネットにつないでドコモお客様サポートを開きます。

「サービス・メールなどの設定」に行きます。

「メール設定」に行けばよさそうに思いますが、ここでハマりました。
spモードを契約した直後は、まだドコモメールが使えないので「各種設定の一覧を見る」に行きます。

ドコモメールはドコモクラウドの機能の一つなので「ドコモクラウド」へ行きます。

「ドコモクラウド設定」へ行きます。

「ドコモメール」が「停止中」になっている場合はそこを押して利用中にします。
途中色々表示されますが画面を保存していません。(それで使えるようになるのか分からず試行錯誤してたので)

ドコモメールが利用中になったら、最初に戻って今度は「メール設定」に行きます。

spモードパスワードを設定していないのに尋ねられたら、初期値は0000です。すぐに変更するように言われます。

「docomo ID利用設定の確認／変更」に行きます。

「docomo IDでドコモメールを利用する」にすることで、spモードアプリ(ドコモ端末にしか入っていない)を使わなくても @docomo.ne.jp のメールアドレスを使えるようになります。

「設定を確定する」これで設定は完了です。

ここまで設定すればドコモメール（ブラウザ版）やIMAPに対応したメールアプリを使ってドコモメールを使えます。
IMAPに対応したメールアプリの設定方法は他のメールソフトからのご利用方法にあります。

spモードを契約すれば他にもdocomo Wi-Fiを使えるようになるのでdocomo Wi-Fiの「0001docomo」をSIMフリー端末で利用し、エリア内に入れば自動接続できるようにするための設定方法。も使えるようになります。
Wi-Fiを使うために毎回Web認証のためのユーザーIDとパスワードを入力する必要がなくなります。

気象データ取得

気象データ取得を完成させて、次回からは取得したデータの時系列解析をします。

区間推定

区間推定について説明しました。

年月を指定して気象データを取得

先週に引き続き、気象庁から取得したデータを、利用可能な形で整形、保存します。

点推定

点推定における推定量の性質について説明しました。

検定の考え方

iPod touchでingress

GPS

スマホ内蔵のGPSより、独立GPSユニットの方が位置の精度が高い筈、と思って試してみました。
Androidの場合、bluetooth接続の独立GPSを使うには「疑似ロケーションを許可」しなければならず、そうするとingressは使えないのでiOSにして、iPhoneは高いし独立GPSユニット使うならGPSを内蔵しないiPod touchでも構わないのでiPod touchにしました。
Appleの認証を受けた独立GPSでなければ、脱獄しないと使えないそうです。認証を受けていて脱獄なしで使えるのはXGPSとBadElfとGNSだけのようです。

iPod touch自身ではWiFiを使った測位しか出来ないので地図を見ると現在位置の誤差を表す円が大きいのですが、独立GPSを繋いで屋外に出ると一気に誤差が小さくなって精度が上がることがわかります。
これでingressやれば細かな位置調整も楽かな、と思ったのですが、そんなに甘くはありませんでした。
精度だけなら普段使っているXperiaZ3よりも高い筈なのですが、それでも建物による反射もあるでしょうし、ポケットに入るような小型GPSの精度はそんなに高くはありません。
スマホの場合は平滑化などの演算で、誤差が大きくてもそれなりに滑らかに移動しているように誤魔化してくれているのでしょう。どうせポータルの位置もずれていることがありますから、本来の位置と少々ずれていても、直前の位置から極端に動かなければingressには問題ありません。
一方、独立GPSの場合、平滑化などの演算を行っていないようで、本来の位置とのずれは小さい筈ですが、そのずれの範囲内で直前の位置とは関係なく唐突に動くため、XperiaZ3と比べてかなりingressを使い辛いです。全く使えないわけではないですが、Xperiaなら簡単に出来る位置調整に、すごく手間取った挙句に出来ないことがあり、Xperaiが壊れるまでiPod touchと独立GPSの組み合わせでingressをすることはないと思いました。
ネットを見ていると、スマホのGPSが良くないので独立GPSを買ったら良くなったという人の話を見ますが、多分スマホの方をちゃんとしたGPSの機種に買い替えたほうが良いです。特にAndroidの場合は機種による差が大きいですから。

タイマーアプリ

Androidの場合は、ingressの画面にタイマーを重ねて表示できるIntegrated Timer For Ingressというとても便利なソフトがありますが、iOSの場合は重ね合わせが出来ません。(恐らく、重ね合わせが出来るようにするとフィッシングなどの悪意を持ったアプリを作りやすくなるからかな？と思います。)
そこで代わりのアプリを探すと、次の3つが見つかりました。
Vidgetsは時間が来た時の通知機能が機能しませんでした。
Portal Timer for Ingressは、タイマーをスタートするためにはウィジェットではなくアプリそのものを切り替える必要があります。
Hacking Timer for Ingressウィジェットとして使うことができるのでアプリを切り替える必要がなく、通知機能もあるのでこれが良いです。

GPSロガーアプリ

ingressとは関係ないですが独立GPSに関連した話で。
独立したGPSは精度は高いですが、GPS衛星からの電波が届かない地下街などでは使えません。
その点、スマホはGPSの電波が届かない地下でも基地局情報を用いて位置を調べることができます。
GPS衛星を使った測位だけでなく、基地局情報を用いた測位もログとして残せるアプリを調べてみたら、GPSロガーAndroidだけが記録できるようです。iOS用ではそのようなアプリは見つかりませんでした。
保存形式はGPX, KMLの他に、GPS衛星測位だけですがmnea形式でも保存できますし、ログを自動的にメールで送ったりdropboxに送る機能もあるのでパソコンへの転送も楽で、GPSのログを記録するアプリとしてはこれが一番良いと思います。
Eye-Fi x Android の組み合わせが超全自動転送&詳細ジオタグ付きにできてすごい便利
kindleFire HDXとEye-Fi mobiカードが作ってくれる最強のフォトストレージ環境

Webクエリ

Excelを用いてインターネットのデータを自動取得する方法を説明しました。

中心極限定理

標本平均の確率分布を正規分布を用いて近似する方法を説明しました。

do while文を使った繰り返し

標本分布

確率変数の独立性と、独立な確率変数の和の期待値、分散について説明しました。

確率分布

統計学において、推定しようとする対象は確率分布です。その確率分布に関する基本的な用語を説明しました。

環境データ解析学

統計学II

ヒストグラムの作り方について説明しました。

スライドがちゃんと表示されない場合は、この記事のタイトルの「統計学II」をクリックしてこの記事だけを表示させてみてください。
講義の後で統計検定について連絡用のメールアドレスを教えてくれた人には、その日の夕方にメールを送りました。
もしも届いていないようでしたら、私宛にメールか直接、連絡ください。

中心極限定理

沢山の観測値の和の確率分布を近似する中心極限定理を説明し、それを用いて実際に母平均を推定する方法を説明しました。

次回の課題は

標本平均と不偏分散

標本平均、標本分散の確率分布について説明しました。

行列の対角化

行列の対角化に必要な固有値、固有ベクトルについて説明しました。

余因子

余因子行列、順列の符号を用いた行列式の公式を説明しました。

確率分布の平均と分散

確率分布の平均と分散について説明しました。

4次以上の行列式

演習

第二章の問題演習をしました。

行列式

行列式の計算方法を説明しました。

正規分布

もっともよく使われる正規分布について説明しました。

主な確率分布

主な確率分布を紹介しました。離散型は二項分布が基本となります。

逆行列

逆行列の求め方を説明しました。

ソルバーを用いた最尤推定量の求め方

今日のテーマ

最尤推定量は(対数)尤度を最大にする値としてデータから計算できます。この計算式θ=θ(X₁,X₂,...,X_n)を具体的に書き表すことが出来る問題もあります(以前のレポートで出しました)が、計算式を求めることが出来ない問題もあります。
そこで今日はExcelを使って最尤推定量を求めて、AICを計算し、モデルを選ぶ方法を勉強しましょう。

準備

Excelを立ち上げたら左上の「ファイル」を押して、「オプション」をクリックします。そして左側の「アドイン」をクリックして、下の「管理:Excelアドイン」の横の「設定」をクリックします。

するとアドインの画面に切り替わるので「ソルバーアドイン」にチェックをつけてOKをクリックします。

サンプルデータ

男児 女児 期間(週) 体重(g) 期間(週) 体重(g) 40 2968 40 3317 38 2795 36 2729 40 3163 40 2935 35 2925 38 2754 36 2625 42 3210 37 2847 39 2817 41 3292 40 3126 40 3473 37 2539 37 2628 36 2412 38 3176 38 2991 40 3421 39 2875 38 2975 40 3231

出典:一般化線形モデル入門, Annette J.Dobson 著, 共立出版
これをコピーしてエクセルに貼り付けてください。散布図を書くとこのようになります。

考えるモデル

パラメータが少ない順番に

男女共通モデル

男女ともY=α+βX+ε, ε～N(0,σ²)
つまり妊娠期間Xに対して体重Yの確率分布は
男女ともN(α+βX,σ²)
推定するパラメータはα, β, σの3つ。

傾きだけ共通モデル

男児:Y=α₁+βX+ε, ε～N(0,σ²)
女児:Y=α₂+βX+ε, ε～N(0,σ²)
つまり妊娠期間Xに対して体重Yの確率分布は
男児:N(α₁+βX,σ²)
女児:N(α₂+βX,σ²)
推定するパラメータはα₁, α₂, β, σの4つ。

男女別モデル

男児:Y=α₁+β₁X+ε, ε～N(0,σ²)
女児:Y=α₂+β₂X+ε, ε～N(0,σ²)
つまり妊娠期間Xに対して体重Yの確率分布は
男児:N(α₁+β₁X,σ²)
女児:N(α₂+β₂X,σ²)
推定するパラメータはα₁, α₂, β₁, β₂, σの5つ。

まず最初に、パラメータが一番少ない男女共通モデルについて、α, β, σをExcelに求めてもらいます。この問題に関しては既に小テストで、α, βをデータを用いて書き表す式を求めていますが、そのような式を求めることが出来ない場合でもExcelに求めてもらうための練習として、この問題に関してもExcelに求めてもらいます。
そのためには、まず最初に適当な初期値を与えます。するとExcelのソルバーはそこから出発して真の値に近づいていきます。初期値をどのように与えればよいでしょうか。最初から真の値が分かるならExcelのソルバーに頼る必要はないのですから、真の値に近くてしかも簡単に求められる値を初期値として使います。
直線y=α+βxのαとβの両方を簡単に知ることは出来ませんが、水平な線、つまりβ=0ならば、その高さαは平均値にしておけば大体データに近そうです。この時σは平均値からのずれ、つまり標準偏差にしておきましょう。
そこで早速体重の平均値と標準偏差を求めましょう。Excelに次のように貼り付けます。

A15に平均、A16に標準偏差と書いてから、B15に平均を求める式、B16に標準偏差を求める式を書きます。
B15にはB3からB14までとD3からD14までの平均ということで =AVERAGE(B3:B14,D3:D14) と書きます。先頭の等号を忘れないでください。
B16にはB3からB14までとD3からD14までの標準偏差ということで =STDEV(B3:B14,D3:D14) と書きます。
するとこのようになります。

セルの幅に合わせて四捨五入されますので、人によって表示される値が少し違うこともあります。
この値を参考にしてα, β, σの初期値を決めます。
まずA17, A18, A19セルにα, β, σと書いて
B17セルにはαの初期値として平均値を書きます。適当なところで四捨五入してかまいません。
B18セルにはβの初期値を書きます。傾きが良く分からなくて水平線を引くので0にします。
B19セルにはσの初期値として標準偏差を書きます。これも適当なところで四捨五入してかまいません。

対数尤度の計算

確率関数や密度関数にデータを代入した値が尤度、それにlogをつけたのが対数尤度です。独立な複数の観測値に関して尤度は掛け算になっているので、logをつけて掛け算を足し算にします。
log(Π_i=1ⁿf(X_i))=Σ_i=1ⁿlog(f(X_i))
正規分布N(μ,σ²)の密度関数は =NORMDIST(x, μ, σ, 0) です。最後の0を1に変えると分布関数になります。
今回の問題では、xのところが体重yで、平均値μがα+βxです。それにlogを付けますが、Excelでは自然対数は =LN(x) ですのでE3セルに =LN(NORMDIST(B3,B$17+B$18*A3,B$19,0)) と書きます。対応するセルに色がつくので確認してください。

α, β, σに相当するB17, B18, B19の数字の前には$を付けています。これはこの式をコピーして下に貼り付けても番号が変わらないようにするためです。
E4セルからE14セルにも同様に数式を書きます。一つ一つ書いていると手間がかかりますのでE3セルをコピーして貼り付けます。数式の中のセルの番号がどのように変わるか確認してください。
次にF3セルに女児の一人目の対数尤度の式 =LN(NORMDIST(D3,B$17+B$18*C3,B$19,0)) を書きます。

これもF4からF14へコピーして、E3からF14まで全部足したものが対数尤度です。
E15セルに対数尤度と書いて、F15セルに合計を求める式 =SUM(E3:F14) を書きます。

ソルバー

ここが今回のポイントです。これからα, β, σに色々な値を代入して対数尤度が最大になればよいのですが、試行錯誤するのは大変ですのでExcelのソルバーを使います。
エクセルの上の列の「データ」をクリックして、右側の「ソルバー」をクリックします。「ソルバー」が見当たらなかったら、このページの上の準備のところを見てください。

するとこのような画面が表示されるので、F15を最大にするためにB17からB19を変化させる、と指定します。
また、αは正の数とは限らないので「制約のない変数を非負数にする(K)」のチェックを外します。さらにソルバーの設定を変えるのでオプションを押します。

「GRG非線形」のタブを選んで、微分係数を「中央」にします。

OKを押してソルバーの画面に戻って「解決(S)」を押すと、対数尤度が大きくなるようにα, β, σの値が調整されて次のメッセージが表示されます。

OKをクリックすると、対数尤度の値が最大になっています。

AICの計算

このようにして求めた対数尤度の最大値からパラメーターの個数を引いた値が一番大きなモデルが、一番良いモデルなのですが、この方法が考えられた経緯から、
-2×(対数尤度の最大値)+2×(パラメータの個数)
をAIC (Akaike information criterion)といい、マイナスを掛けているのでAICが一番小さなモデルが一番良いモデルです。
このモデルのAICは-2×(-159.265)+2×3=324.53です。

レポート課題

傾きだけ共通モデル、及び男女別モデルについて、最尤推定量及びAICの値を求め、どのモデルが一番良いか選びなさい。