「日々の学び」一覧

同じ出処から取ってきたデータを全て訓練データとして使わずに、 訓練データとテストデータに分割して、訓練データで作ったモデルに対するテストデータの精度を返す、 みたいなことをやるらしい。交差検証(CrossValidation)という名前が付いている、 sklearn.model_selection.cross_val_score( estimator, X, y=None, groups=None, scoring=None, cv=’warn’, n_jobs=None, verbose=0, fit_params=None, pre_dispatch=‘2*n_jobs’, error_score=’raise-deprecating’) estimatorとして、例えば決定木分類であればDecisionTreeClassifierのインスタンスを渡す。 Xは説明変数、yは目的変数。交差検証自体のアルゴリズムを選択できてcvに渡す値で制御できる。 学習済みのモデルを渡してスコアが戻るのではなく、 単なるモデルのインスタンスと学習前のデータを別々に渡している作りなのを見ると、 モデル毎のスコアを並べてみたくなる..。そういう使い方するのか? \"スコア\"と言ってるこの値は具体的には何の値なのか..? K-分割交差検証(K-fold cross-validation)の場合cvは以下のようになる。 cvを省略するとK=3が使われる。使ったデータは例のあやめ。 clf = DecisionTreeClassifier(max_depth = 3) # integer, to specify the number of folds in a (Stratified)KFold, score1 = cross_val_score(estimator = clf, X = iris.data, y = iris.target, cv = 10) # None, to use the default 3-fold cross validation, score2 = cross_val_score(estimator = clf, X = iris.data, y = iris.target) 決定木の深さを交差検証で求める(失敗) suumoから引いてきた家賃~占有面積,築年数,階数,バストイレ別データについて、 目的変数を家賃、説明変数を占有面積、築年数、階数、バストイレ別として決定木を作ってみた。 決定木の深さは交差検証で求める、と書かれているので以下のようにしてみた。 import sys import pandas as pd import matplotlib.pyplot as plt from sklearn import tree from sklearn.model_selection import train_test_split from sklearn.model_selection import cross_val_score path = \"./fuchu.csv\" train = pd.read_csv(filepath_or_buffer=path) train = train.drop_duplicates() feature_name = [\"area\",\"age\",\"bt_separated\"] feature_name2 = [\"area\"] train_x = train[feature_name2].values train_y = train[\"value\"].values sort_idx = train_x.flatten().argsort() max_depth_list = [1,3,5,7,9,11,13,15,17] for max_depth in max_depth_list: regressor = tree.DecisionTreeRegressor(max_depth = max_depth) regressor.fit(train_x,train_y) score = cross_val_score(estimator=regressor,X=train_x,y=train_y) print([max_depth,score.mean()]) 結果、おかしい...。マイナスってなんだ。 [1, -0.020456885057637583] [3, 0.012376893473610817] [5, -0.014519973352621932] [7, -0.03332646095737912] [9, -0.06862181989491711] [11, -0.09590667631145984] [13, -0.1345799052512994] [15, -0.1529266053016934] [17, -0.17054978124438266] 説明変数を占有面積だけに絞って散布図書いて、その上に回帰結果をプロットすると、 異常値の0に引っ張られて乱れてた...。以下は深さ=17の時。むぅ. 真ん中の塊の中を左下から右上に向かってギザギザに進んで欲しいのだが、 物凄い引っ張られよう。異常値が原因というよりは、引っ張られ過ぎなんだな。

[mathjax] アンサンブル学習とかランダムフォレストに入門する前に決定木に入門する。 決定木はやっていることが直感的でわかりやすい。 決定木回帰と決定木分類。 ここよりはドメインとの連結部分が大変なんだろうと思った。 あと、Pythonは練習しないとな。 CART(Classification and Regression Tree)法(単に分類と回帰を英語にしただけだ！) 木を作るのだけれども、それが面白かったので今回と次回で書いてみる。 決定木回帰 (y=(x-0.5)^2)という2次曲線に従う事象があるとして、(y)を観測するとする。 観測による誤差が平均(mu=0)、分散(sigma^2=0.1)の正規分布に従うとして(y=(x-0.5)^2+epsilon)。 区間([0.0,1.0])の間に等間隔に存在する観測値(x,y)。 試行毎に(epsilon)が変わってくるので、毎回異なる。 この区間に(16)個のデータがあるとして、それを訓練データとして使ってモデルを作る。 (x_{train},y_{train})とする。 モデルの作成（学習）はfit()。 Scikit-learnに全て用意されていてデータを放り込むだけでモデルが出来る。 決定木の葉の最大値を(5)としている。他のパラメタは全部デフォルト値。 import numpy as np # 区間[0,1]上に16個の点を等間隔に生成する X_train = np.linspace(start=0,stop=1,num=16) y_train = (X_train - 0.5) ** 2 + np.random.normal(loc=0.0,scale=0.1,size=16) # 16x1配列を1x16に整形 X_train = X_train.reshape(16,1) print(X_train) # 決定木回帰 from sklearn.tree import DecisionTreeRegressor DTR = DecisionTreeRegressor(max_leaf_nodes=5) DTR.fit(X_train,y_train) 出来上がったモデルにテスト用データを流し込んでみる。 区間([0.0,1.0])に100個のデータを発生させてpredict()を呼ぶ。 最後、訓練データと回帰結果を同じグラフを書いてみて終了。 # 区間[0,1]上に100個の点を等間隔に生成する X_test = np.linspace(0,1,100) X_test = X_test.reshape(100,1) # 回帰! y_predict = DTR.predict(X_test) X_train = X_train.reshape(16) X_test = X_test.reshape(100) # 描画 import matplotlib.pyplot as plt plt.scatter(X_train,y_train) plt.plot(X_test,y_predict) plt.savefig(\'img.png\') plt.show() 葉の最大値を(5)としたので、木の深さが規定されて、 階段の個数が決まっている。 以下、6回分のモデルと予測の同時プロット。 (epsilon)の変化により訓練データが微妙に変わるだけで、 決定木の構造がむちゃくちゃ変化するのが特徴。 それっぽく言うとロバスト性が無いとか。 訓練データによって決定木がかなり違うことを利用して、 複数の決定木から多数決で結果を得ようというのがアンサンブル学習の試み。 むー。順々に詰めていこう。 max_leaf_nodesをデータの個数-1と一緒にすると、 全ての訓練データを通るモデルを作ることができる。 訓練データに対しては100%の精度が出るが、未学習のデータに対して答えられなくなる(ほんとに？)。 これが過学習(overfitting)。 Rによる 統計的学習入門posted with amazlet at 19.04.22Gareth James Daniela Witten Trevor Hastie Robert Tibshirani 朝倉書店 売り上げランキング: 118,792Amazon.co.jpで詳細を見る

[mathjax] 今回から、atomにmathjax-wrapperを入れてatomで数式入り文章を書いてみる。 数式を書くと瞬時にプレビューができて、10倍は速くなった気がする。1時間くらいで書いてみる。 母集団の(sigma)、(mu)が既知なのなら、共役事前分布は正規分布で良いのだから、 一体何をやりたいのか...という記事になってることに1日経って気づいた...。 チラシの裏程度なので気にしない。 何故、事前分布と事後分布を同じ確率分布に揃えるかというと、計算済みの事後分布を次の事前分布として 使うことができるから。事後分布の計算を繰り返していくと理論上は精度が上がっていく。 事前分布と事後分布を同じ分布にするために、逆ガンマ分布を事前分布として採用する。 そうすることで、左辺の事後分布と、右辺の尤度x共役事前分布が共に逆ガンマ分布となる。 (母平均が既知で分散が未知、という条件がつく。) 正規分布の共役事前確率分布として逆ガンマ分布を使う、というところまでは良かったのだけども、 逆ガンマ分布のパラメタが巷には2種類なのと4種類なのがあってもやもや。 ガンマ分布をさらに一般化した一般化ガンマ分布が4パラメータなのだそうな。 そもそもガンマ分布はカイ二乗分布と指数分布の一般化なので、どこまで一般化するのか。 ひとまず触れてはいけないところに触れてしまったようなので、2パラメータで出直してみる。 合ってるんだか間違ってるんだかも不明なのだけども、結論としては、細かいことはどうでもよくて、 事後分布と尤度x共役事前分布の関係を脳裏に焼き付けるためのプロセス。 何周かしないと真理にはたどり着けない...。 ガンマ分布と逆ガンマ分布 ガンマ分布の確率密度関数は以下。(alpha)は形状パラメータ,(beta)はスケールパラメータ。 (alpha) を固定して (beta)を動かすとカイ二乗分布。 (beta) を固定して (alpha)を動かすと指数分布。 Excelで　(alpha) と (beta)をグリグリ動かしてみると意味がわかる。面白い。 良くできてるなー。 $$ begin{eqnarray} f(x|alpha,beta) = frac{1}{beta^alpha Gamma(alpha)} x^{alpha-1} expleft(-frac{x}{beta}right) end{eqnarray} $$ そして、(alpha) を大きくしていくと形状が正規分布に近づいていく。 累積確率分布関数は以下。 $$ begin{eqnarray} F(x|alpha,beta) = frac{1}{beta^alphaGamma(alpha)} int_0^xt^{alpha-1} exp (e^{-frac{t}{beta}})dt end{eqnarray} $$ 対して逆ガンマ分布の確率密度関数は以下。ガンマ分布で(x\'=frac{1}{x}) とすると出てくる。 $$ begin{eqnarray} f(x|alpha,beta)=frac{beta^{alpha}}{ Gamma(alpha)} x^{-(alpha+1)} expleft(-frac{beta}{x}right) end{eqnarray} $$ The shorthand for the distribution, X~inverted gamma(α,β), or IG(α, β), means that a random variable X has this distribution with positive parameters α and β. 再びベイズ統計へ ベイズの定理を使って出てくる事後確率分布と事前確率分布の関係は以下の通りだった。 $$ begin{eqnarray} p(mu,sigma^2|boldsymbol{x})&=&frac{p(boldsymbol{x}|mu,sigma^2)p(mu,sigma^2)}{p(boldsymbol{x})} \\ &propto& p(boldsymbol{x}|mu,sigma^2)p(mu,sigma^2) end{eqnarray} $$ 尤度は周辺確率の積として表されるけれども、母集団の平均が既知で分散が未知の場合の話をするので、 確率変数(X)を分散にもつ正規分布を考える。(y)は標本平均。(mu)は既知の母平均。 $$ begin{eqnarray} p(x | mu, X) &=& N(x|mu,X) \\ &=& frac{1}{sqrt{2pi x}}expleft(-frac{(y-mu)^2}{2x}right) end{eqnarray} $$ 正規分布である尤度に逆ガンマ分布をかけて式変形していく。 比例の関係を見たいので定数をのぞいてシンプルにするところがポイント。 $$ begin{eqnarray} p(mu,sigma^2|boldsymbol{x}) &propto& frac{1}{sqrt{2pi x}}expleft(-frac{(y-mu)^2}{2x}right) frac{beta_0^alpha}{Gamma(alpha_0)} left(frac{1}{x}right)^{alpha_0+1} expleft(frac{-beta_0}{x}right) \\ &propto & frac{1}{x} exp left( - frac{(y-mu)^2}{2x}right) x^{-alpha_0-2} exp left( -frac{beta_0}{x} right) \\ &propto & x^{-alpha_0-2} exp left( -frac{1}{x}(beta_0 + frac{1}{2}(y-mu)^2 right) end{eqnarray} $$ 右辺の確率分布が逆ガンマ分布であるためには、定数のぞいた箇所、つまり(x)の肩と(exp)の肩が 逆ガンマ分布のそれと同じでなければならないから、 $$ begin{eqnarray} -(alpha_0-2) &=& alpha +1 \\ alpha &=& alpha_0+frac{1}{2} \\ end{eqnarray} $$ また、 $$ begin{eqnarray} beta &=& beta_0 + frac{1}{2} (y-mu)^2 end{eqnarray} $$ これで、左辺の逆ガンマ分布の (alpha)、(beta)が作れる。 更新するたびに (alpha)、(beta)が増加していくのがわかる。

[mathjax] 2項分布に従う母集団の母平均を推測するために有意水準を設定して95%信頼区間を求めてみた。 母平均のあたりがついていない状況だとやりにくい。 [clink url=\"https://ikuty.com/2019/01/11/sampling/\"] (hat{p})がどんな値であっても下限は(hat{p})の関数で抑えられると思ったので、 気になって(hat{p})を変数のまま残すとどうなるかやってみた。 begin{eqnarray} 1.96sqrt{frac{hat{p}(1-hat{p})}{n}} le 0.05 \\ frac{1.96}{0.05}sqrt{hat{p}(1-hat{p})} le sqrt{n} \\ 39.2^2 hat{p}(1-hat{p}) le n end{eqnarray} 左辺を(f(hat{p}))と置くと (f(hat{p}))は下に凸の2次関数であって、 (frac{d}{dhat{p}}f(hat{p})=0)の時に最大となる。というか(hat{p}=0.5)。 (hat{p}=0.5)であるとすると、これはアンケートを取るときのサンプル数を求める式と同じで、 非常に有名な以下の定数が出てくる。 begin{eqnarray} 1537 * 0.5 (1-0.5) le n \\ 384 le n end{eqnarray} (hat{p})がどんな値であっても、サンプル数を400とれば、 有意水準=5%の95%信頼区間を得られる。 だから、アンケートの(n)数はだいたい400で、となる。 さらに、有意水準を10%にとれば、(n)の下限は100で抑えられる。 なるはやのアンケートなら100、ちゃんとやるには400、というやつがこれ。

[mathjax] ちょっと不思議な計算をしてみる。 仮定に仮定を積み重ねた素人の統計。 成功か失敗かを応答する認証装置があったとする。 1回の試行における成功確率(p)は試行によらず一定でありベルヌーイ試行である。 (n)回の独立な試行を繰り返したとき、成功数(k)を確率変数とする離散確率変数に従う。 二項分布の確率密度関数は以下の通り。 begin{eqnarray} P(X=k)= {}_n C_k p^k (1-p)^{n-k} end{eqnarray} 期待値、分散は、 begin{eqnarray} E(X) &=& np \\ V(X) &=& np(1-p) end{eqnarray} (z)得点（偏差値,つまり平均からの誤差が標準偏差何個分か？）は、 begin{eqnarray} z &=& frac{X-E(X)}{sigma} \\ &=& frac{X-E(X)}{sqrt{V(X)}} \\ &=& frac{X-np}{sqrt{np(1-p)}} end{eqnarray} であり、(z)は標準正規分布に従う。 これを標本比率(hat{p}=frac{X}{n})を使うように式変形する。 begin{eqnarray} z &=& frac{frac{1}{n}}{frac{1}{n}} frac{X-np}{sqrt{np(1-p)}} \\ &=& frac{frac{X}{n}-p}{sqrt{frac{p(1-p)}{n}}} \\ &=& frac{hat{p}-p}{sqrt{frac{p(1-p)}{n}}} end{eqnarray} (n)が十分に大きいとき、(z)は標準正規分布(N(0,1))に従う。 従って、(Z)の95%信頼区間は以下である。 begin{eqnarray} -1.96 le Z le 1.96 end{eqnarray} なので、 begin{eqnarray} -1.96 le frac{hat{p}-p}{sqrt{frac{p(1-p)}{n}}} le 1.96 end{eqnarray} (hat{p})は(p)の一致推定量であるから、(n)が大なるとき(hat{p}=p)とすることができる。 begin{eqnarray} -1.96 le frac{hat{p}-p}{sqrt{frac{hat{p}(1-hat{p})}{n}}} le 1.96 \\ end{eqnarray} (p)について解くと(p)の95%信頼区間が求まる。 begin{eqnarray} hat{p}-1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} le p le hat{p}+1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} end{eqnarray} 上記のにおいて、標準誤差(1.96sqrt{frac{hat{p}(1-hat{p})}{n}})が小さければ小さいほど、 95%信頼区間の幅が狭くなる。この幅が5%以内であることを言うためには以下である必要がある。 (有意水準=5%) begin{eqnarray} 1.96sqrt{frac{hat{p}(1-hat{p})}{n}} le 0.05 end{eqnarray} 観測された(hat{p})が(0.9)であったとして(n)について解くと、 begin{eqnarray} 1.96sqrt{frac{0.9(1-0.9)}{n}} le 0.05 \\ frac{1.96}{0.05} sqrt{0.09} le sqrt{n} \\ 11.76 le sqrt{n} \\ 138.2 le n end{eqnarray} 139回試行すれば、100回中95回は(p)は以下の95%信頼区間に収まる。 つまり95%信頼区間は以下となる。 begin{eqnarray} hat{p}-1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} &le& p le hat{p}+1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} \\ 0.9-1.96 frac{sqrt{0.09}}{sqrt{139}} &le& p le 0.9 + 1.96 frac{sqrt{0.09}}{sqrt{139}} \\ 0.9-1.96 frac{0.3}{11.78} &le& p le 0.9+1.96 frac{0.3}{11.78} \\ 0.85 &le& p le 0.95 end{eqnarray} (n)を下げたい場合は有意水準を下げれば良い。 統計的に有意水準=10%まで許容されることがある。 有意水準が10%であるとすると、(n)は35以上であれば良いことになる。 begin{eqnarray} 1.96sqrt{frac{0.9(1-0.9)}{n}} le 0.1 \\ frac{1.96}{0.1} sqrt{0.09} le sqrt{n} \\ 5.88 le sqrt{n} \\ 34.6 le n end{eqnarray} 信頼区間と有意水準の式において(p)を標本から取ってきたけど、 アンケートにおいてYes/Noを答える場合、(p)は標本における最大値（つまり0.5)を 設定して(n)を求める。 つまり、(p)として利用するのは標本比率ではないのかな？と。 このあたり、(hat{p})を変数として残すとどういうことがわかった。 [clink url=\"https://ikuty.com/2019/01/13/sampling_with2/\"]

[mathjax] いつか起こる大地震がもし昨日起きたとしたら明日起こる確率は下がるのか？ 飛行機が昨日落ちたとしたらしばらくは飛行機は落ちないのか？ うまくいかない人生が今日またうまくいかなかったとして将来うまくいかない確率は下がるのか？ 起こるか起こらないかの確率が変わらないのであれば、将来は過去に影響されないらしい。 影響されるかどうか、と聞かれるとされない、と答えるだろうけど、それを説明することができるらしい。 そんな無記憶性の件読んでみたのでまとめてみる。 統計学入門とこちらを参考にさせて頂きました。 幾何分布 幾何分布 二項分布、ポアソン分布は(n)回の試行のうち(x)回事象(A)が発生したときを話題にしていた。 前提を少し変えて、予め試行する回数を決めないで、 (x)回目の試行で初めて事象(A)が発生した、 という話をすることもできる様子。 確率変数(x)が等間隔に並ぶ時刻であるとすることで、 事象(A)が発生するまでの待ち時間に関する確率分布を作れる。 (x)回目の試行で初めて事象(A)が起こった、ということを、 (x-1)回事象(bar{A})が起こり、次に事象(A)が起きたと考える。 事象(A)の生起確率を(p)、事象(bar{A})の生起確率を(q)とすると、 その確率は以下のようになる。 begin{eqnarray} f(x) = p cdot q^{x-1} end{eqnarray} (x)の増加1回に対して(q)を1回かける構造で、公比(q)の等比数列になっている。 等比数列って英語でgeometric seriesって言うもんだから、幾何分布っていう名前がついてる様子。 (f(x))の確率変数(x)は事象(A)が発生するまでの試行回数（時間）である、と読む様子。 では確率変数(X)が幾何分布に従うとき、期待値はどうかというと(frac{1}{p})となる。 確率の以下の読み方から、(x)回の試行に平均して(frac{1}{p})回かかる、というのは妥当に見える。 事象(A)の生起確率(p)の読み方について、その逆数(frac{1}{p})は、事象(A)が起こるまでの試行回数と読む。 (p)回の試行で初めて事象(A)が起こった、というシーンで事象(A)の生起確率を(p)と置いているので..。 期待値 幾何分布の期待値は(frac{1}{p})。ベルヌーイ試行の確率が(p)であるならば、平均(frac{1}{p})回で事象が起こる。 分散は(frac{1-p}{p^2})。期待値の証明は以下の通り。へぇ。 begin{eqnarray} E(X) &=& sum_x x cdot f(x) \\ &=& sum_x x cdot p q^{x-1} end{eqnarray} (E(X)がfrac{1}{p})であることを示したい。 恒等式を使ってやる奴ではなく、愚直にやる奴を書く。 まず、(q=1-p)として(E(X))を変形しておいてスタート。 begin{eqnarray} E(X) &=& sum_x x cdot p (1-p)^{x-1} \\ &=& p sum_x x cdot (1-p)^{x-1} end{eqnarray} 右辺を生み出すために、(frac{1}{x-1})のテイラー展開を持ち出す。 begin{eqnarray} frac{1}{1-x} &=& 1 + x + x^2 + cdots \\ &=& sum_{k=0}^{infty} x^k end{eqnarray} 左辺を(x)で微分すると以下の通り。 begin{eqnarray} left( frac{1}{1-x} right) frac{d}{dx} = frac{1}{(1-x)^2} end{eqnarray} 右辺を(x)で微分すると以下の通り。 begin{eqnarray} sum_{k=0}^{infty} x^k frac{d}{dx} = sum_{k=1}^{infty} k x^{k-1} end{eqnarray} なので、 begin{eqnarray} frac{1}{(1-x)^2} = sum_{k=1}^{infty} k x^{k-1} end{eqnarray} (x=1-p)として式変形すると、 begin{eqnarray} frac{1}{p^2} = sum_{k=1}^{infty} k (1-p)^{k-1} end{eqnarray} (E(X))にこれらを代入すると、 begin{eqnarray} E(X) &=& p sum_x x cdot (1-p)^{x-1} \\ &=& p cdot frac{1}{p^2} \\ &=& frac{1}{p} end{eqnarray} 本当に(frac{1}{p})になった。 両辺を微分したものが等しいって、なんでだっけ？ 無記憶性 どうも、世の中には(n-1)回連続して失敗して(n)回目で初めて成功することを言っているものと、 (n)回連続して失敗して、(n+1)回目で初めて成功することを言っているものがある。 期待値も分散も若干違うものになる。 ここからは(n)回の失敗に続いて(n+1)回目で初めて成功するケースに切り替える。 その時の確率を(P(X=n))とする。で、失敗が(n)回以上連続して起こる確率(P(Xgeq n))を考える。 begin{eqnarray} P(X geq n) &=& P(X=n) + P(X=n+1) + P(X=n+2) + cdots \\ &=& p(1-p)^n + p(1-p)^{n+1} + p(1-p)^{n+2} + cdots \\ &=& p(1-p) left( 1 + (1-p)^1 + (1-p)^2 + cdots right) \\ &=& p(1-p) sum_{k=1}^{infty} (1-p)^{k-1} end{eqnarray} 途中の無限級数は(frac{1}{1-x})の級数展開になっていて、 以下みたいになる。 begin{eqnarray} P(X geq n) &=& p(1-p)^n sum_{k=1}^{infty} (1-p)^{k-1} \\ &=& p(1-p)^n frac{1}{1-(1-p)} \\ &=& p(1-p)^n frac{1}{p} \\ &=& (1-p)^n end{eqnarray} (n)回連続して失敗した上で、さらに連続して(k)回の失敗を重ねる確率を考える。 (n)回連続して失敗する確率は(P(Xgeq n))。 この条件の上でさらに(k)回失敗を重ねる確率は条件付き確率として(P(Xgeq n+k | X geq n))。 条件付き確率の定義と乗法定理から式を展開していく。（ここが難しかった...) begin{eqnarray} P(Xgeq n+k | X geq n) &=& frac{P((X geq n+k)cap (X geq n) )}{P(X geq n)} \\ &=& frac{P(X geq n+k)}{P(X geq n)} \\ &=& frac{(1-p)^{n+k}}{(1-p)^n} \\ &=& (1-p)^k \\ &=& P(X=k) end{eqnarray} ということで、以下が成り立つことがわかる。 begin{eqnarray} P(Xgeq n+k | X geq n) &=& P(X=k) end{eqnarray} よーく見てみると、(n)回連続して失敗した後に(k)回連続して失敗する確率と、 (n)回の失敗無しに、最初から(k)回連続して失敗する確率が同じである、と言っている。 凄まじいことに、(n)回連続して失敗することは、次の(k)回の失敗に全く影響を及ぼさない、と言っている。 何回失敗しようと次に失敗する確率はこれまでの失敗に影響されない。 つまり失敗する確率は過去の影響を受けない。 美しすぎる感じがする。

[mathjax] 期待値、分散、歪度、尖度...、確率分布を形成する確率密度関数の特徴を表す値で、 実は、相互に変換できる値なのだという...。読んでいったら若干感動したのでまとめてみる。 40近いオッさんがはじめてテイラー展開のありがたさを味わう瞬間の記録。 モーメント母関数とモーメント モーメント母関数を以下のように定義。 begin{eqnarray} M_x(t) &=& E(e^{tX}) \\ &=& int_{-infty}^{infty} e^{tx} f(x) dx end{eqnarray} 英語で書くとmoment generating function。モーメントを作る関数。 ここで(f(x))というのが確率密度関数。 もともとこの時点で積分が存在しないかもしれない。 確率密度関数によっては、期待値、分散、歪度、尖度を直接求めるのは難しい。 しかし、モーメント母関数の(r)階導関数からモーメント(mu_r)を解析的に求められる性質から、 期待値、分散、歪度、尖度を求めることができる。 さて...、40近いオッさんは思い出す。 指数関数のテイラー展開は、マクローリン級数を使って以下の通り。 begin{eqnarray} e^x = 1+x+x^2/2!+x^3/3!+cdots end{eqnarray} (tX)だと、 begin{eqnarray} e^{tX} = 1+tX+(tX)^2/2! +(tX)^3/3! + cdots end{eqnarray} 両辺の期待値をとる。右辺全体の期待値はそれぞれの項の期待値の和にできるので、 begin{eqnarray} E(e^{tX}) &=& M_x(t)\\ &=& E(1) + E(tX) + E((tX)^2/2!) + E((tX)^3/3!)) + cdots end{eqnarray} (t)に対する定数を出すと begin{eqnarray} M_x(t) &=& E(e^{tX})\\ &=& E(1) + E(X)t + E(X^2/2!)t^2 + E(X^3/3!)t^3 + cdots \\ &=& 1 + mu_1 t + (mu_2/2!)t^2 + (mu_3/3!)t^3 + cdots end{eqnarray} キター。おわかりだろうか...。 (M_x(t))は(t)に関する展開式の係数に各次数のモーメントを含んでいる。 (t)について1回微分すると0次までの項は消える。 2次以降の項の(t)の次数が1減って残る。1次の項だけ(t)が消える。 そのとき(t=0)とすると、階数の係数(M^{(r)}_X(0)=mu_r)だけ残る！ つまり、以下のようなとんでもないことになる。 begin{eqnarray} M_X\'(0) = mu_1 \\ M_X\'\'(0) = mu_2 \\ M_X\'\'\'(0) = mu_3 \\ end{eqnarray} 各次数のモーメントである期待値、分散、歪度、尖度を、 モーメント母関数の(r)階導関数から求められるということになる。 指数分布のモーメント 試しに指数分布でやってみる。 begin{eqnarray} M_x(t) &=& int_{0}^{infty} e^{tx} lambda e^{-lambda x} dx \\ &=& lambda int_{0}^{infty} e^{(t-lambda)x} dx end{eqnarray} 指数関数の積分のところでおっ、と思ったけど、以下となる。 begin{eqnarray} M_x(t) = frac{ lambda }{ lambda -t} end{eqnarray} これ、解析的に微分できるのかな...と思うんだけども高校数学で暗記するやつ。 微分と積分を行ったり来たりするとわかる。 begin{eqnarray} M_x^{(1)}(t) &=& frac{ lambda }{ (lambda -t)^2} \\ M_x^{(2)}(t) &=& frac{ 2 cdot lambda }{(lambda -t)^3} \\ M_x^{(3)}(t) &=& frac{ 2 cdot 3 cdot lambda}{(lambda -t)^4} end{eqnarray} (t=0)とおくと、 begin{eqnarray} mu_1 &=& frac{1}{lambda} \\ mu_2 &=& frac{2}{lambda^2} \\ mu_3 &=& frac{2 cdot 3}{ lambda^3} \\ mu_4 &=& frac{2 cdot 3 cdot 4}{ lambda^4} end{eqnarray} これで、微分の数値計算をしなくても解析的に(mu_1)から(mu_4)が求まった。 そして永遠に微分し続けることで指数分布を形作る指標が決まっていく。 すごいなぁ...。