「日々の学び」一覧

[mathjax] 線形単回帰で推定する回帰係数の精度を評価する方法を読んだのでまとめてみる。 当然、真の直線はわからないのだけれども、真の直線があると仮定した上で 推定した回帰係数との関係を考えることで、回帰係数の精度について話せるようになる。 回帰係数の導出 データポイントが(n)個ある状況。 ( (x_1,y_1),(x_2,y_2),cdots,(x_n,y_n) ) 回帰係数(hat{beta_0})と(hat{beta_1})を使って線形回帰したい。 begin{eqnarray} hat{y} = hat{beta_0} + hat{beta_1} x end{eqnarray} データポイントと回帰直線の差を残差平方和(RSS,redisual sum of square)で表す。 データポイントは既に与えられているデータなので、(hat{beta_0},hat{beta_1})の関数。 begin{eqnarray} f(hat{beta_0},hat{beta_1}) = (y_1 -hat{beta_0}-hat{beta_1}x_1)^2 + (y_2 - hat{beta_0}-hat{beta_1}x_2)^2 + cdots + (y_n - hat{beta_0}-hat{beta_1}x_n)^2 end{eqnarray} RSSを最小にする(hat{beta_0})と(hat{beta_1})を求めるために、(hat{beta_0})、(hat{beta_1})それぞれで偏微分して(0)として解く。 なんでそれぞれ個別に偏微分して0と置いて良いかは、 RPML読もうとして力尽きたときに理解したので省略。 参考にした本に( hat{beta_0}),(hat{beta_1}),RSSの3次元の図があって、確かにそれで良さそうな予感。 begin{eqnarray} frac{partial}{partial hat{beta_0}} f(hat{beta_0},hat{beta_1}) = 0 \\ frac{partial}{partial hat{beta_1}} f(hat{beta_0},hat{beta_1}) = 0 \\ end{eqnarray} 以下のようになるらしい。(bar{x})、(bar{y})はデータポイントの標本平均。 なので、データポイントがわかれば計算で求まる。 begin{eqnarray} hat{beta_1} &=& frac{sum_{i=1}^n (x_i-bar{x}) (y_i-bar{y}) }{sum_{i=1}^n (x_i-bar{x})^2 }\\ hat{beta_0} &=& bar{y}-hat{beta_1}bar{x} end{eqnarray} 母回帰直線の推定 データポイントが同じであれば(hat{beta_0}),(hat{beta_1})は同じになるけれども、 データポイントを取り直して異なるデータセットにすると、(hat{beta_0}),(hat{beta_1})は微妙に違う値になる。 じゃあ、データセットを大量に用意したとして、(hat{beta_0}),(hat{beta_1})を計算しまくると、 どこかに収束するんじゃなかろうか。 標本が大量にあると標本平均は母平均に収束する。標準偏差はより小さくなる。 つまりデータが大量にあると、母平均からのズレが小さくなっていく。 大数の弱法則、中心極限定理、ルートnの法則。 begin{eqnarray} hat{sigma} &=& frac{sigma}{sqrt{n}} \\ hat{sigma}^2 &=& frac{sigma^2}{n} end{eqnarray} begin{eqnarray} lim_{n rightarrow infty} hat{sigma}^2 = lim_{n rightarrow infty} frac{sigma^2}{n} = 0 end{eqnarray} [clink url=\"https://ikuty.com/2018/07/17/sample_sigma/\"] (hat{beta_0}),(hat{beta_1})は母回帰直線からどれくらいばらついているのか。 (hat{beta_0}),(hat{beta_1})の分散は以下を使うらしい。 両方に出てくる(sigma^2)は、母回帰直線と回帰直線の差となる項の散らばり度合い。 つまり、(Y=beta_0 + beta_1 X + epsilon )としたときの(epsilon)の分散。 begin{eqnarray} sigma_{hat{beta_0}}^2 &=& sigma^2 Bigl[frac{1}{n} + frac{bar{x}^2}{sum_{i=1}^n (x_i-bar{x})^2} Bigr] \\ sigma_{hat{beta_1}}^2 &=& frac{sigma^2}{sum_{i=1}^n (x_i -bar{x})^2} end{eqnarray} (x_i)が散らばれば散らばるほど、(sigma_{hat{beta_1}}^2)は小さくなる。 データポイントの(x)成分が小さい方から大きい方まで含まれれば、傾き(beta_1)を推定しやすくなる。 そして、(bar{x}=0)であるならば、(hat{beta_0})の散らばりは、(hat{mu})の散らばりと等しくなる。 最終的に求めたいのは不明な(sigma^2)だが、(sigma^2)はデータから計算できる。 (sigma)の推定値(RSE,Resual Standard Error)はRSSから推定する。 begin{eqnarray} sqrt{frac{f(hat{beta_0},hat{beta_1})}{(n-2)}} end{eqnarray} (hat{beta_1})の標準偏差がわかったので、95%信頼区間を求めることができる。 線形回帰における(hat{beta_1})の95%信頼区間は、 begin{eqnarray} Bigl[ hat{beta_1} - 1.96 sigma_{hat{beta_1}},hat{beta_1} + 1.96 sigma_{hat{beta_1}} Bigr] end{eqnarray} 同様に(hat{beta_0})の95%信頼区間は、 begin{eqnarray} Bigl[ hat{beta_0} - 1.96 sigma_{hat{beta_0}},hat{beta_0} + 1.96 sigma_{hat{beta_0}} Bigr] end{eqnarray}

稼働中のEC2を落とさないでALB下で切り替えた作業記録を書いてみます。 こちら↓の方が詳しく書いてあります..。今回書いた記事の特徴は以下となります。 AutoScalingグループを使わない ALBの下で切り替える deploy手順をAWSの外に用意する [clink implicit=\"false\" url=\"https://qiita.com/keitakn/items/6abe6c971e4dec3b69ef\" imgurl=\"https://camo.qiitausercontent.com/08bed869c98443e0474ce8ce78bdbe964a09f1e9/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f37313839392f62303437313164612d623035362d313262662d646164612d3439653930333231663366342e6a706567\" title=\"AWS CodeDeploy でEC2のBlue/Greenデプロイを作成する\" excerpt=\"AWS CodeDeploy を使ってBlue/Green デプロイの仕組みを構築する為の手順を紹介します。Blue/Greenデプロイとは？現在稼働している環境と別にもう1つ稼働環境を作成し、ロードバランサー等のルーティングを新環境に向けるデプロイ方法です。常にリクエストを受けている稼働中のサーバを置き換えるよりも安全にデプロイ可能なのがメリットになります。\"] [arst_adsense slotnumber=\"1\"] 現状と動機 WPCoreアプデとセキュリティパッチはマメに当てないといけないと痛感して、 仕方なくBlue Green的な方法を導入してみた。 ALBの下にWebサーバ1台。Webサーバの下にDB用EC2が1台(非RDS)。 アップデート対象はWebサーバのみ。 Webサーバ内でWordPressが12個、SNIで動いてる。 x.smallクラス。CloudWatchを見るとLoadAverageは常時30%くらい。 全てのサイトで、pluginでファイルとDBをS3にバックアップしている。 localで開発/WPCore,plugin update/動作確認後、止めずにansibleでdeploy。 deploy実行前に\"メンテナンス中\"に設定。deploy完了後に解除する。 アップデート中は管理画面操作禁止を通達。 deploy、パッチ当てでコケると、S3から戻すまで止まる! S3から戻らないと終わる。 やったこと deploy、アップデート時にのみ、WebサーバのEC2をコピーする。 ALBのターゲットグループにコピーしたEC2を追加する。 元のEC2に対してansibleでdeployする。 元のEC2に対してパッチアップデートする。 ALBの先を元のEC2に戻す。 コピーしたEC2を削除する。 メディアライブラリにファイルをアップロードすると差分が発生するため、 元EC2と一時EC2のファイル達を同期させないといけないけれど、 メンテ中は、管理画面操作を禁止できるという状況であることと、 もともとEC2が1台なのでその仕組みを作っていないから、それは2台以上に増えたときに..。 AMI作成 まずAMIを作成。AMIとはAmazon Machine Imageの頭文字。 [clink implicit=\"false\" url=\"https://docs.aws.amazon.com/ja_jp/AWSEC2/latest/UserGuide/ec2-instances-and-amis.html\" imgurl=\"https://docs.aws.amazon.com/ja_jp/AWSEC2/latest/UserGuide/images/architecture_ami_instance.png\" title=\"インスタンスと AMI\" excerpt=\"Amazon マシンイメージ (AMI) は、ソフトウェア構成 (オペレーティングシステム、アプリケーションサーバー、アプリケーションなど) を記録したテンプレートです。AMI から、クラウドで仮想サーバーとして実行される AMI のコピーであるインスタンスを起動します。以下の図に示すように、1 つの AMI の複数のインスタンスを起動することができます。\"] 手順は以下。 ダッシュボードからコピーしたいインスタンスを選択 アクション->イメージ->イメージの作成を選択 デフォルトだと、AMI作成時にコピー元インスタンスが自動的に停止し、コピー後に自動的に再起動する。 \"再起動しない\"にチェックをいれることで、コピー元インスタンスの停止/再起動を抑止できる。 \"再起動しない\"にチェックをいれないでAMIを作成すると、コピー元が止まってしまうので注意! イメージの作成を押下すると、作成処理が始まる。 EBSが30GiBだと、完了まで1時間程度要してしまった。 ダッシュボード -> AMI から AMIの作成状況を確認できる。 ステータスが available になれば完了。 インスタンスの起動 作成済みAMIからインスタンスを起動する。 ダッシュボード -> AMI を開く 起動したいAMIを選択する アクション -> 起動を押下 すると、インスタンスタイプを聞かれる。 進捗状況はEC2ダッシュボードで確認できる。 ALBのターゲットグループ変更 既にALBのターゲットグループに元EC2が属していて、 セキュリティグループが正しく設定済みで、ヘルスチェックが通っている前提。 現状、ALBの下は元EC2だけなので、AvailabilityZoneは1種類だけ。 ダッシュボード -> ターゲットグループ そこに、新しく作成したインスタンスを追加する。 新しいインスタンスのセキュリティグループを旧インスタンスに合わせて ALBからのInboundを受けられるようにすること。 新しいインスタンスのヘルスチェックが無事通って2台構成になった図。 元のEC2をターゲットグループから削除 元のEC2をターゲットグループから削除する。 EC2ダッシュボードのモニタリングタブをみて、CPU使用率などに変動があることを確認する。 元のEC2に対してゴニョゴニョする ALBのターゲットグループから外れた元のEC2に対してdeployなりパッチ当てを実行する。 元のEC2にElasticIPを当てておけば、再起動してもIPアドレスは変わらない。 この手順においては、新規作成したEC2にElasticIPを当てる必要はない。 元のEC2のセキュリティグループを0.0.0.0/0:80からアクセスできるようにして、 hostsにElasticIPを書いてアクセスするなどの方法で、元のEC2にアクセスする。 このためだけにALBを作って置いておくとそのALBに対して課金されてしまう。 新規作成したEC2を無料枠でやったりしてスペックが低いとパフォーマンスが下がる。 非稼働のEC2に当てたElasticIPで課金される。 一時EC2をターゲットグループから外し、元のEC2を投入する 上の手順を逆から実施。つまり、一時EC2をターゲットグループから外し、元のEC2を投入する。 [arst_adsense slotnumber=\"1\"]

[mathjax] 損失関数をバイアス項、バリアンス項、削減不能誤差の和に分解できることと、 損失は削減不能誤差より下回らないこと、バイアス項、バアリアンス項のトレードオフが起こること、 を読んだ。過学習っていうのはこういうことなのか、と腑に落ちたので記念に記事を書いてみる。 （式変形は細かいところで間違ってるのと、おっさんのチラシの裏なので参考にしないでください） 2乗損失の期待値の式変形 モデルを作った後、訓練データ、テストデータそれぞれの全データについて、 2乗損失の期待値(MSE)を求め、モデルの当てはまりの良さを調べるらしい。 2乗損失を以下のように式変形する。条件付き期待値(E(t|x))ってなんだ...。 begin{eqnarray} L(y(x),t)^2 &=& (y(x)-t)^2 \\ &=& (y(x)-E(t|x)+E(t|x)-t)^2 \\ &=& left( left( y(x)-E(t|x)right) + left( E(t|x) - y(x)right) right)^2 \\ &=& (y(x)-E(t|x))^2 + 2(y(x)-E(t|x))(E(t|x)-t) + (E(t|x)-t)^2 \\ end{eqnarray} 2乗損失の期待値(MSE)は以下。 第2項は(x)、(t)で積分するとゼロになる! begin{eqnarray} E[L(y(x),t)^2] &=& E[ (y(x)-E(t|x))^2 + 2(y(x)-E(t|x))(E(t|x)-t) + (E(t|x)-t)^2 ] \\ &=& E[ (y(x)-E(t|x))^2 + (E(t|x)-t)^2 ] end{eqnarray} 和の期待値は期待値の和なので、 begin{eqnarray} E[L(y(x),t)^2] = E[ (y(x)-E(t|x))^2 ] + E[(E(t|x)-t)^2 ] end{eqnarray} (x)の出処がテストデータではなく訓練データですよ、と明示するために、 以下みたいな書き方に改める。この式の中で(y(x;D))が学習で得られるモデル。 第2項は学習とは関係なく発生する数値。 begin{eqnarray} E_D[L(y(x;D),t)^2] &=& E_D[ (y(x;D)-E(t|x;D))^2 ] + \\ && E_D[(E(t|x;D)-t)^2 ]] end{eqnarray} 第1項の式変形を続ける。括弧が多すぎて力尽きた..。 余計な項を足して引いて次の式変形の足しにするタイプ。 begin{eqnarray} E_D[ (y(x;D)-E(t|x;D))^2 ] &=& E_D[ ( { y(x;D)-E_D[(y(x;D))] } ] &+& { E_D[ y(x;D)] - E[t|x;D])^2 } \\ &=& E_D [ { (y(x;D))-E_D[ y(x;D)] }^2 ] + \\ &=& E_D [ { E_D[ y(x;D)-E[t|x;D] ] }^2 ] end{eqnarray} バイアス・バリアンスと削減不能誤差 以下はバリアンス項と書かれている。 モデル((y(x;D))による予測が訓練データ集合によって変動する度合いの期待値。 異なる訓練データを使ったときにどの程度モデルが変化するかを表す。 過学習の度合い。 begin{eqnarray} E_D bigl[ bigl{ (y(x;D))-E_D[ y(x;D)] bigr}^2 bigr] end{eqnarray} 以下はバイアス項と書かれている。 複雑な事象を単純なモデルで近似したことによる誤差、と書かれてる。 例えば、3次関数+ノイズから発生するデータを直線で近似すると、モデルが単純すぎて値が大きくなる。 モデルが複雑になればなるほどバイアス項は減っていく様子。 未学習の度合い。 begin{eqnarray} E_D bigl[ bigl{ E_D[ y(x;D)]-E[t|x;D] bigr}^2 bigr] end{eqnarray} で、一番最初に出てきたモデルと関係ない以下。 バイアス、バリアンス共に非負の値だから、2乗損失の期待値は以下より小さくなることはない。 奇跡的にバイアス、バリアンス共にゼロだったとしても、以下は学習とは関係なく発生する。 削減できない誤差。 begin{eqnarray} E_Dbigl[bigl(E(t|x;D)-tbigr)^2 ]bigr] end{eqnarray} 結局よくわからない...。体感の結論.. 訓練データを使ってモデルを複雑にしていけばいくほど、 モデルが訓練データにフィットするようになるが、 その訓練データにフィットしまくったモデルは、未知のテストデータを予測しづらくなる。 モデルの複雑度が\"ある程度\"のところまでは、バリアンスの上昇よりもバイアスの低下が効くから、 訓練データに対する2乗誤差、テストデータに対する2乗誤差ともに減少する。 モデルの複雑度が\"ある程度\"を超えると、バイアスの低下が頭打ちになる一方でバリアンスが上昇し、 訓練データに対する2乗誤差が低下する一方で、テストデータに対する2乗誤差が上昇する。 どう頑張っても、削減不可能な誤差が存在する。 条件付き期待値(E(t|x))の意味を理解できずプロットすることは叶わなかった。

[mathjax] おっさんが入門した軌跡シリーズです。損失関数に関して学んだことをメモしておきます。 入力値(x)、正解(t)からなる訓練データ(T=(x,t))が大量に与えられたときに、 (f(x,w))によって回帰なり分類なりをする。 仮に立てた(f(x,w))と正解(t)の距離(Lleft(f(x,w),tright))を損失関数と呼んだり。 あえて(Lleft(f(x,w),tright))としているのは、 一番わかりやすそうな残差の2乗だけでなく、他があるから。 2乗損失 残差の2乗だと、(f(x,w))と(t)の差が大きい場合に必要以上に大きくなってしまう。 ほとんどのデータで残差が0なのに、特殊なデータで残差が100とかになられたら、 全体の損失は測れそうにないし、異常値に敏感すぎる。 begin{eqnarray} Lleft(f(x,w),tright) = (t-f(x,w))^2 end{eqnarray} Huber損失 単に残差の2乗を使うだけでは不十分で、\"(f(x,w))と(t)の差の大小にあまり影響されないこと\"が必要。 残差の絶対値がある値を超えるまでは残差の2乗、 超えてからは線形にするというのもある(Huber損失)。 begin{eqnarray} Lleft(f(x,w),tright) = begin{cases} (t-f(x,w))^2 & t in [t-delta,t+delta] \\ 2delta (|t-f|-frac{delta}{2}) & それ以外 end{cases} end{eqnarray} 損失関数が微分可能か 損失関数を最小(極小)にすることが目的なので..、 損失関数の1階導関数(勾配ベクトル)を使ってパラメタを更新したりする。 begin{eqnarray} w_{i+1} = w_i - eta left(f\'(x,w),tright) end{eqnarray} 損失関数が微分可能（連続）だと勾配ベクトルがすぐに求まるので、 損失関数は不連続点をなくすように作るらしい。 (Huber損失の2乗から線形に切り替わるところは連続になってる。) 次回以降実際のデータとモデルを使ってやってみる。

[mathjax] 集合に対して再帰的に境界を入れていく操作が決定木の作成。 では、集合のどこに境界を入れれば良いか。 属性をテストすることにより得られる情報量が最も大きくなるように入れる。 汎化能力、みたいな言葉を読んでいくにあたってこの先、結構抽象的な話になるので一度確認。 データが(n)個のクラスに分類できるとして、クラス(i)に属する確率を(P_i)とする。 このとき、あるデータがクラス(i)に属することを知るには(log_2 frac{1}{P_i})の情報量が必要。 その期待値は(I(P_1,P_2,cdots,P_n)=sum_{i=1}^{n} P_i log_2 frac{1}{P_i} = - sum_{n}^{i=1} P_i log_2 P_i)。 情報量の平均をエントロピーとか。 (D_p)個のデータを(D_{left})、(D_{right})に分割するとする。 その時、属性(f)に関する問いを使って分割する。 2分木の左ノード、右ノードだからleft,right。 (I(D_p))は分割前のエントロピー。 (I(D_{left}))、(I(D_{right}))は分割後のエントロピー。 分割前には(I(D_p))ビット必要だったのが、 分割後には(frac{N_{left}}{N_p} I(D_{left}) + frac{N_{right}}{N_p} I(D_{right}))ビット必要になった。と読むらしい。 その差を情報ゲイン(IG(D_p))と呼んで以下のように定義するらしい。 begin{eqnarray} IG(D_p,f) = I(D_P) - frac{N_{left}}{N_p} I(D_{left}) - frac{N_{right}}{N_p} I(D_{right}) end{eqnarray} 分割前よりも分割後の方が(IG(D_p,f))だけエントロピーが低い、という事実に関して、 分割に使った問いにより、(IG(D_p,f))の情報量を獲得した、と考えるらしい。 情報ゲインが最大になるような問いを根にもってきて、 再帰的に情報ゲインが大きいものから順に問うことで決定木を作っていく。

同じ出処から取ってきたデータを全て訓練データとして使わずに、 訓練データとテストデータに分割して、訓練データで作ったモデルに対するテストデータの精度を返す、 みたいなことをやるらしい。交差検証(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/\"]