2020-01-02

ポアソン回帰とモデル選択

R データサイエンスモデル選択モデリング GLM

新年明けましておめでとうございます．
前回G検定対策の記事を書いて以降，論文を書いていたりなので全然更新ができていませんでした．G検定は無事合格したので，今後は更新頻度を上げて，E資格のための記事も書いていこうと思います．

今回は下記の本の3章について取り上げたいと思います．

データ解析のための統計モデリング入門一般化線形モデル・階層ベイズモデル・MCMC （シリーズ確率と情報の科学） [ 久保拓弥 ]

ジャンル: 本・雑誌・コミック > 科学・医学・技術 > 数学
ショップ: 楽天ブックス
価格: 4,180円

統計モデリング入門

植物100個体の種子数のデータを用います. データの詳細として下記に６個体のデータを示します.
また，
y：種子の数
x：個体の大きさ,
f：肥料を与えた場合T，肥料を与えてない場合にC
です．

> d<-read.csv("data.csv")

   y     x f
1  6  8.31 C
2  6  9.44 C
3  6  9.50 C
4 12  9.07 C
5 10 10.16 C
6  4  8.32 C

個体の大きさや肥料の影響は種子数の数に影響あるのかな??

またどのような統計モデルを当てはめれば，データの特徴をより表すことができているかについて考えます．
f:id:koki12070930:20191118223755p:plain

個体の大きさが影響する場合

データの性質上，ポアソン分布を仮定します．
ある個体iの平均種子数を $\lambda_i$ とします．
ポアソン分布のパラメータ $\lambda_i$ は
$\lambda_i=\mathrm{exp}(\beta_1+\beta_2x_i)$
となります．
ただし， $x_i$ は個体 $i$ の大きさです．
$x_i$ の項があることから, 個体の大きさを考慮してることになっています. パラメータ $\beta_2$ の値が大きければ, 平均種子数は大きくなるように $\beta_1, \beta_2$ の値はとても重要であることが分かります. これらパラメータ $\beta_1, \beta_2$ を求める代表的な方法は最尤法です. Rではこの計算を自動でやってくれるパッケージがあるので今回はそれを用います!

> fit1<-glm(y~x,data=d,family=poisson(link="log"))
> fit1

Call:  glm(formula = y ~ x, family = poisson(link = "log"), data = d)

Coefficients:
(Intercept)            x  
    1.29172      0.07566  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    89.51 
Residual Deviance: 84.99 	AIC: 474.8

したがって,
$\lambda_i=\mathrm{exp}(1.29172+0.07566 x_i)$ と求めることができます.
このパラメータを用いてポアソン分布をを可視化します．
f:id:koki12070930:20200102001122p:plain

施肥処理が影響する場合

因子型説明変数はダミー変数に変換します．
\begin{eqnarray}
d_i&=&
\begin{cases}
0~(f_i=C)\\
1~(f_i=T)
\end{cases}
\end{eqnarray}
$\lambda_i=\mathrm{exp}(\beta_1+\beta_3d_i)$ であるので， $\beta_1$ ， $\beta_3$ のパラメータ推定をします．

> fit2<-glm(y~f,data=d,family=poisson)
> fit2

Call:  glm(formula = y ~ f, family = poisson, data = d)

Coefficients:
(Intercept)           fT  
    2.05156      0.01277  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    89.51 
Residual Deviance: 89.48 	AIC: 479.3

したがって,
肥料なし $\to \lambda_i=\mathrm{exp}(2.05156+0)$
肥料あり $\to \lambda_i=\mathrm{exp}(2.05156+ 0.01277 \times 1)$ 　

個体の大きさと施肥処理が影響する場合

ここでは $\lambda_i=\mathrm{exp}(\beta_1+\beta_2x_i+\beta_3d_i)$ です．

>fit3<-glm(y~x+f,data=d,family=poisson)
> fit3

Call:  glm(formula = y ~ x + f, family = poisson, data = d)

Coefficients:
(Intercept)            x           fT  
    1.26311      0.08007     -0.03200  

Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
Null Deviance:	    89.51 
Residual Deviance: 84.81 	AIC: 476.6

したがって,
肥料なし $\to \lambda_i=\mathrm{exp}(1.26311+0.08007x_i)$
肥料あり $\to \lambda_i=\mathrm{exp}(1.26311+0.08007x_i -0.032)$ 　
となります．

先ほど同様，これらを用いてポアソン分布を可視化します．
f:id:koki12070930:20200102003328p:plain

当てはまりが良いモデルはどれ??

ここまで，様々な状況下でパラメータ推定を行なって来ました．
結局どのモデルが一番当てはまりが良いのでしょうか??

最尤法は，対数尤度の最大化問題を考えてパラメータ推定を行うのでした．
では問題です．

これまでのモデルの中で対数尤度が最大のものが一番良いモデルであると結論づけて良いでしょうか??

答えは"良くない"です．
モデルの選択の基準として，最大対数尤度を用いることは適切ではなく，AICを用いたモデル選択が重要になります．
なぜ最大対数では適さないのかは下記の記事を参考にしてください．
koki12070930.hatenablog.com

今回はモデルの選択には情報量基準AICを推奨しましたが，目的によっては尤度比検定を行うことがあります．
ポアソン分布以外のGLMの話題や混合モデルへの拡張について，次回以降に書きたいと思います!

ではまた!

2019-09-23

ディープラーニングG検定対策

機械学習 G検定 AI ディープラーニング人工知能

G検定に出てくであろう重要な用語をピックアップしたいと思います．

人工知能とは
人工知能の歴史
人工知能の研究の問題

人工知能とは

人工知能(Artificial Intelligence)とは，コンピュータによる学習，推論，など人間の知能の働きを実現するための研究分野．ただ「何が人工知能なのか」研究者や専門家の中でも共有されてた明確な定義はありません.

機械学習(Machine Learning)とは，「人工知能」を実現するための手法の一つです．
アーサーサミュエルによる「明示的にプログラムしなくても学習する能力をコンピュータに与える研究分野」というのが有名です．

深層学習(Deep Learning)とは，機械学習アルゴリズムの中の一つです．ニューラルネットワークを多層にしたモデルになっています．

まとめると，人工知能>機械学習>深層学習です！！

人工知能の歴史

人工知能の研究はこれまで，3度のブームを迎えてきました．一つずつ解説します．

・第一次AIブーム (推論探索の時代)
簡単な迷路やパズルのような問題(トイ・プロブレム)を解くことに成功しました．これでも当時では活気的なことでしたが，複雑な問題は解くことができず冬の時代を迎えることになります．

・第二次AIブーム (知識の時代)
データベースに「知識」を大量に放り込んだ「エキスパートシステム」と呼ばれるシステムがこの当時に活躍しました．しかし，知識を蓄積，管理することが困難であることが明らかになりまたもや冬の時代を迎えることになりました．

・第三次AIブーム (機械学習と深層学習の時代)
人工知能自らが知識を獲得する機械学習が実用化されました．つまり，第二次AIブームの問題点を克服したのです．さらに特徴量を人工知能が自ら獲得する深層学習が台頭したことが第三次AIブームのきっかけと言えます．

人工知能の研究の問題

フレーム問題

1969年に提唱された人工知能の研究の最難解である問題です．
「今しようとしていることに関係のある事柄だけを選び出すことが，実は非常に難しい」という問題です．分かりやすく解釈すると次のようになります．

本来は枠(フレーム)を定め，関係のないことを排除する必要があります．
しかし人工知能は，関係のあること(フレーム内)と関係のないこと(フレーム外)を判断することが困難です．

フレーム問題に縛られないAI＝強いAI
フレーム問題に縛られたAI=弱いAI

シンボルグラウンディング問題

シンボルグラウンディング問題とは．記号とその対象の意味がいかにして結びつくかという問題です．
我々は，「火」の意味も「トカゲ」の意味も理解しているのでポケモンで出てくる「ヒトカゲ」を初めて見たとしても，「あいつが噂のヒトカゲでは？」と推測できます．しかしコンピュータはヒトカゲは記号の羅列にすぎないので「火のついたトカゲ」と記述できたとしてもその意味を理解することができません．

このように記号とその意味するものが結びつかないことをシンボルグラウンディング問題と言います．
ちなみにヒトカゲのくだりは僕の大好きなレペゼン地球の動画の中で取り上げられています．勉強になりました．リンクを貼っておきます．ついでにチャンネル登録よろしくです．

【ポケモンの名前なんて見ただけで当てれるわw】天才が挑む...!!!『DJ社長レペゼン地球』

シンギュラリティ(技術的特異点)

シンギュラリティとは「人工知能が人間の知能を超えて人間に変わって文明の主役にとってかわる」という問題です．

レイ・カーツワイルは「シンギュラリティは2045年に到来する」と述べています．
超越的な人工知能がもたらす影響に対し，警鐘を鳴らす人々も数多くいます．マイクロソフトの創業者であるビルゲイツは「私も人工知能に懸念を抱く側にいる1人だ」と脅威論に賛同しています．
これら脅威論に対し日本では，2014年に人工知能学会に倫理委員会が設置されました．

G検定に出てくるであろう覚えるべき細々した用語を少しまとめました．
参考書として以下がとても役に立ちます．

深層学習教科書ディープラーニング G検定（ジェネラリスト）公式テキスト（EXAMPRESS） [ 一般社団法人日本ディープラーニング協会 ]

価格: 3080 円
楽天で詳細を見る

徹底攻略ディープラーニングG検定ジェネラリスト問題集 [ 明松真司 ]

価格: 2310 円
楽天で詳細を見る

11月のG検定はさっさと合格して，2月のE試験に備えたいと思います．ブログの登録よろしくお願いします．
ではまた！

2019-09-02

FiFAランキングのデータから分かるブラジルの凄さ

Python データ分析ライブラリ数学 pandas データサイエンス soccer

本記事では，kaggleで提供されているFIFA Soccer Ranking*1のデータセットを用いた簡単な考察を行います．

CSVファイルの読み込み
日本のFIFAランキングの推移
ブラジルの凄さ

CSVファイルの読み込み

まずはcsvファイルの読み込みです．

#データの読み込み
import pandas as pd
df=pd.read_csv("fifa_ranking.csv")

次にデータの先頭５行を表示します．

df.head() 

   rank country_full     ...     confederation   rank_date
0     1      Germany     ...              UEFA  1993-08-08
1     2        Italy     ...              UEFA  1993-08-08
2     3  Switzerland     ...              UEFA  1993-08-08
3     4       Sweden     ...              UEFA  1993-08-08
4     5    Argentina     ...          CONMEBOL  1993-08-08

[5 rows x 16 columns]

1993年8月では，1位がドイツ，2位がイタリア，3位がなんとスイス!だったみたいです．

次に日本のFIFAランキングに注目したいと思います．

data=df[df["country_full"]=="Japan"]
rank=data["rank"]
rank.describe() 

count    286.000000
mean      35.048951
std       13.066002
min        9.000000
25%       24.000000
50%       35.000000
75%       45.750000
max       62.000000

日本のFIFAランキングの平均はおおよそ35位であり，1番いい時が9位???

data[data["rank"]==9] 

      rank country_full     ...     confederation   rank_date
8128     9        Japan     ...               AFC  1998-02-18
8320     9        Japan     ...               AFC  1998-03-18

1998年2~3月では，日本のFIFAランキングは9位だったようです．こんなにランキングが良かったカラクリは当時のFIFAランキングのシステムにあるようです．詳しくはこの記事*2
を見てください．

日本のFIFAランキングの推移

次に日本のFIFAランキングの推移の様子を可視化したいと思います．

import matplotlib.pyplot as plt 
data["rank_date"]=pd.to_datetime(data["rank_date"])   
date=data["rank_date"]                 
fig,ax=plt.subplots()          
ax.plot(date,rank) 
plt.show()

f:id:koki12070930:20190902175734p:plain — 日本のFIFAランキングの推移

ブラジルの凄さ

ブラジルについても見ていきたいと思います．
めんどくさいので可視化のためのコードは関数化します．

def fifatrans(country):
    df["rank_date"]=pd.to_datetime(df["rank_date"])     
    data=df[df["country_full"]==country]
    rank=data["rank"]   
    date=data["rank_date"]   
    fig,ax=plt.subplots()          
    ax.plot(date,rank)  
    plt.show() 

fifatrans("Brazil")

f:id:koki12070930:20190902180414p:plain — ブラジルのFIFAランキングの推移

ブラジルやば...ずっと1位だったんや．
箱ひげ図も書いてみます．

def ranking(country):
    data=df[df["country_full"]==country]
    rank=data["rank"]
    return rank

rank1=ranking("Brazil")
rank2=ranking("Germany")
rank3=ranking("Argentina")

fig,ax=plt.subplots()         
ax.boxplot([rank1,rank2,rank3], labels=['Brazil', 'Germany', 'Argentina'])  
plt.show()

f:id:koki12070930:20190902183600p:plain
ドイツやアルゼンチンと比較してもブラジルの凄さが分かります．
今日は疲れたのでここまでにしておきます．明日以降再度編集するかもしれません!

下記の本はnumpyやpandas 等についてかなり丁寧に書いてあるので基礎から学び方にはおすすめです!

Pythonによるあたらしいデータ分析の教科書／寺田学／辻真吾／鈴木たかのり【1000円以上送料無料】

ジャンル: 本・雑誌・コミック > PC・システム開発 > プログラミング > その他
ショップ: bookfan 2号店楽天市場店
価格: 2,678円

ではまた!

*1:https://www.kaggle.com/tadhgfitzgerald/fifa-international-soccer-mens-ranking-1993now

*2:https://www.footballchannel.jp/2013/10/11/post9521/

2019-06-18

ランダム行列の無線通信技術への応用

MIMO 多変量解析固有値ランダム行列

自己情報量
平均情報量
条件付き自己情報量
相互情報量
通信路容量
MIMOチャネルへの拡張

　情報量は, 頻繁に起こる事柄は情報としてあまり価値がなく，滅多に起こらない事柄に関しては情報として価値が高いはずです．情報量をを定量的に表したものとして次の自己情報量があります．

自己情報量

ある事象 $a$ が生起する確率を $P$ とする．このとき，事象 $a$ が生起したことで得られる情報量 $I(p)$ を(自己)情報量と呼び以下の式で定義します．
\begin{eqnarray*}
I(a)=-\mathrm{log}~P(a)
\end{eqnarray*}
$\mathrm{log}$ の底を2にした場合は情報量の単位をbitと呼びます．

平均情報量

n個の排反事象 $a_1,\dots,a_n$ からなる事象系Xを考えます．この時，自己情報量の期待値として，平均情報量を次のように定義します．
\begin{eqnarray*}
\mathcal{H}(X)=\sum_{i=1}^{n}P(a_i)I(a_i)
\end{eqnarray*}
この自己情報量の平均値をエントロピーとも呼びます．

条件付き自己情報量

　以下では，通信回線を通して情報を伝送する場合について議論します．
また，送信シンボルを $x_i$ とする事象系を $X$ ，受信シンボルを $y_i$ とする事象系 $Y$ について考えます．このとき， $y_i$ が受信された条件のもとで正しい送信シンボル $x_i$ を知ることによって得られる条件付き自己情報量 $I(x_i\backslash y_i$ )は次のように定義します．
\begin{eqnarray*}
I(x_i\backslash y_i）=-\mathrm{log}~P(x_i\backslash y_i)
\end{eqnarray*}

相互情報量

相互情報量 $I(x_i;y_i)$ は次のように定義されます．
\begin{eqnarray*}
I(x_i;y_i)=I(x_i)-I(x_i\backslash y_i）
\end{eqnarray*}
相互情報量が意味することは，自己情報量-条件付き自己情報量の差です．
自己情報量と同様に，平均情報量に対して条件付き平均情報量 $\mathcal{H}(X\backslash Y$ )を定義することが可能です．

平均相互情報量は，次のように定義されます．
\begin{eqnarray*}
\mathcal{H}(X;Y)=\mathcal{H}(X)-\mathcal{H}(X\backslash Y)
\end{eqnarray*}
平均相互情報量はシンボルあたりに伝送される情報量の平均値を表しています．つまり一秒間あたりで考えれば，通信路で伝送された情報の伝送速度に対応しています．

通信路容量

　これまでは離散型の場合について議論してきましたが，連続型の情報量の場合に関しては，これまでの議論を $\sum$ を $\int$ に変更すればOKです．下記のものは連続型の場合に呼び名が変わります．

平均情報量 $\to$ 差分エントロピー

条件付き平均情報量 $\to$ 条件付き差分エントロピー

Shannonは単位時間あたりに伝送できる最大の情報量，すなわち，先ほど定義した平均相互情報量の最大値をその通信路(チャネル)の容量と定義しました．

平均相互情報量は
\begin{eqnarray*}
\mathcal{H}(X;Y)&=&\mathcal{H}(X)-\mathcal{H}(X\backslash Y)\\
&=&\mathcal{H}(Y)-\mathcal{H}(Y\backslash X)
\end{eqnarray*}
と与えられましたが，ここで我々がコントロールすることができるのは，送信シンボルの確率密度関数 $p(x)$ のみであることに注意します．
したがって，チャネルの容量は
\begin{eqnarray*}
C=\underset{p(x)}{\mathrm{max}}~\mathcal{H}(X;Y)
\end{eqnarray*}
で与えられます．

MIMOチャネルへの拡張

　携帯電話や無線LANなどの技術が広く社会に普及するに伴い，次世代の通信システムとして無線通信の技術のMIMOに近年期待が高まっています．MIMOでは送信側，受信側で複数のアンテナを利用情報を伝達します．MIMOチャネル(通信路)は電波が不規則な散乱や反射を受けるため，MIMOの通信路モデルを表現するにはランダム行列が重要になります[図1]．

f:id:koki12070930:20190615224638p:plain — 図1

MIMOの数学モデルは下記のように表すことができます．
\begin{eqnarray*}
\bf{y}=\bf{Hx}+\bf{v}
\end{eqnarray*}

MIMOチャネルの容量を複素ウィシャート行列 $\bf{HH'}$ の固有値 $\lambda_i$ を用いて導出できることをTalatar (1999)で示されました．
電力制約 $P$ ，送信シンボル数 $n$ ，受信シンボル数 $m$ に対して，MIMOチャネルの容量は
\begin{eqnarray*}
E\bigg[\sum_{I=1}^{n}\mathrm{log}~\mathrm{det}(1+(P/n)\lambda_i)\bigg]
\end{eqnarray*}
と表されます．

ウィシャート行列の固有値分布を用いて，このチャネル容量を求めることが理論的には可能ですが，固有値分布は数値計算が困難なためモンテカルロ法を用いたシミュレーションの方法が主流となっています．

以上のようにランダム行列は統計学の枠を超えた応用例を持っていることが分かります．
　最近情報幾何学や学習理論の方面でもランダム行列が使われていることを耳にするので，そちらの方についても今後調べていきたいと思います．

今回参考になった本を紹介します．

わかりやすいMIMOシステム技術

作者:大鐘武雄,小川恭孝
出版社/メーカー: オーム社
発売日: 2009/06/01
メディア: 単行本

こちらの本はTelatarの論文の解説や情報理論の基本事項が丁寧に解説されていてとても参考になりました．

ランダム行列の数理と科学

作者:渡辺澄夫,永尾太郎,樺島祥介,田中利幸,中島伸一
出版社/メーカー: 森北出版
発売日: 2014/04/17
メディア: 単行本（ソフトカバー）

こちらの本にもランダム行列の応用例としてMIMOの話題が取り上げられています．

Telatar, I. E. (1999), Capacity of multi-antenna Gaussian channels

ではまた！

2019-06-08

モデルの評価と選択(1)

Python 機械学習データサイエンスモデル選択クロスバリデーション多変量解析

　昨年，某大学の授業を履修していたのですが，最近そこで勉強したことがかなり役に立ち，復習を兼ねてまとめることにします．

どの説明変数をモデルの中に取り入れたら，現象を有効に予測し，モデルを構築することができるでしょうか?この問に，答えるためにはモデル評価基準について学ぶ必要があります．モデル評価基準は，モデルの良さを測る尺度であり，下記にあるのが有名どころかなと思います．

情報量基準AIC (Akaike 1973，1974)
ベイズ型モデル評価基準 BIC(Schwarz，1978)
クロス・バリデーション (Stone，1974)

今回はクロス・バリデーションの必要性について述べたいと思います．

多項式回帰モデルなどにおいては，次数の決定が本質的であり，次数選択が重要なりますが，適切な次数を選択するための注意点について述べたいと思います．

多項式回帰モデルと残差平方和

$x_i:$ 説明変数， $\hat{y}_i:$ 予測値, $\hat{\beta_i}:$ 回帰係数
予測値(推定した多項式モデル)　
$\hat{y}_i=\hat{\beta}_0+\hat{\beta}_1x_i+\hat{\beta}_2x_i^2+\cdots+\hat{\beta}_px_i^p$ ，
for $i=1,\cdots,n$
これに対して，残差 $e_i$ は，観測データ $y_i-$ 予測値 $\hat{y}_i$ すなわち， $e_i=y_i-\hat{y}_i$ で求められ，
残差の二乗和(残差平方和)を ${\displaystyle \begin{eqnarray} \sum_{i}^{n}e_i^2=(y_i-\hat{y}_i)^2 \end{eqnarray} }$ と表します．
ここで問題です．

多項式モデルの次数を上げていくと，残差平方和はどのようになるでしょうか?
線形回帰モデルの場合に，説明変数の個数を増やしていくと残差平方和はどのようになるでしょうか?

$\rightarrow$ 単調に減少し，次第に0に近くが正解です．

残差平方和はモデルのデータへの適合度の良さを測っているように見えるが，最適なモデルを決定する基準としては機能していません．では最適な次数を決めるにはどのようにすればよいか？
構築したモデルの良さは予測の観点から評価しなければなりません．
モデルを推定するために用いた学習データとモデルを評価するために用いるテストデータは分ける必要があります．残差平方和がモデルを決定する基準として適切ではないことを下記の図で示します．

f:id:koki12070930:20190608193122p:plain — 多項式の次数

上記の図は， $y=\sin(2\pi x)$ に標準偏差0.3を加え，そこからサンプリングしたデータを用い，最小二乗法で得られた多項式と $y=\sin(2\pi x)$ との残差平方和の推移を表しています．

過学習

　図から分かるように，次数が上がるにつれて，学習データに基づいたモデルと学習データとの残差平方和の値は小さくなります．しかし，学習データに基づいたモデルとテストデータとの差つまり予測の観点からの残差平方和は次数が4の時が最適で，それ以降は誤差が増加する傾向にあります．テストデータに対する「予測力」は，次数が4以上では高まらないことを示しています．このように学習データのみにチューニングした状況を過学習(オーバフィッティング)と言います．

構築したモデルの良さは予測の観点から評価するべきであることが分かりました．しかし統計的モデリングの過程において，観測測定されたデータに加えて，別個に将来のデータを獲得するのは困難です．今回のように，学習データとテストデータに分けてモデルを評価すればいいのですが，テストデータを増やせば，そのぶん学習データを減らすことになります．
この問題を解決するのがクロス・バリデーションです．次回はこのクロス・バリデーションについて説明します．

今回は下記の本に書いてることにも一部取り上げました．Pythonで実際に動かしながら理解できるので大変勉強になります．