2015-12-15

自転車ポエム

このエントリはpyspa Advent Calendar 2015の15日目の記事です。pyspa Advent Calendarは「自分が好きなものや興味があるものをなんでもいいので書く」という事で、数年に一度あるか無いかのポエムです。

ロードバイクを入手してからの1年

2014年冬にCannondale SYNAPSE CARBON 5を購入、弱虫ペダルで言うと手嶋先輩が乗ってる奴。そして2015年は輪行を修得した事で行動範囲が一気に拡大した。輪行とは専用の袋に自転車をバラして格納して鉄道などの交通機関に乗る技である。
出張先での練習が可能に

2015年で印象に残った道

荒川サイクリングロード

自宅から5分でアクセスできたのもあり、一番良く走ったのがここ。葛西臨海公園と河口から40km付近の秋ケ瀬公園の間が普段の練習コース。
実際にはサイクリングロードでは無く緊急用河川敷道路として整備されているため舗装が良く幅員も広い。都内最強の平坦道ではなかろうか。左岸を堀切橋あたりから河口まで行くと片道10kmだが、途中にゲートが無く、ランナーとサイクリストしかいないためペースキープの練習にもってこい。日没後は闇に紛れた全身黒ずくめのランナーが出没するため危険。

江戸川自転車道

首都圏外郭放水路で有名な龍Q館に行くにはこの道が便利。江戸川の始点となる利根川も含めて都心まで100km走って関東平野の広大さを実感した。難点としては単調な道が続き景色に変化が無いため、10分で飽きる。いたる所で拡張工事がされているため、今後の整備に期待できそう。

「海まで127.5kmです」の看板、500m毎に立っていて親切?

渡良瀬川

利根川から入り群馬の市街地を経由して遡上して足尾銅山まで行った。観光地なので登りの練習をしながらも道の駅に寄ったり、渋い店構えの蕎麦屋で補給できる。わたらせ渓谷も美しい。足尾に差しかかった辺りからSIRENの舞台を彷彿とさせる昭和な世界が広がっており、味わいがあった。

あらゆる鉄が錆びている

三浦海岸

三浦海岸から三崎口へ抜けて城ケ崎公園へ。砂浜、キャベツ畑、風力発電機、魚屋と変化に富んだ牧歌的なコースで和む。幸運にも風が穏かな日に走れた。オススメのおみやげは鮪のみりん干し。

都民の森

その柔らかな響きとは裏腹に、斜度10%越えのつづら折りに殺されかけた難所。ヒルクライムの苦しさと共に楽しさに目覚めた場所でもある。

白馬 - 木崎湖周辺

北アルプス山麓グランフォンドの120kmコースに出場。生憎の雨でDNF(Don’t Finish)だったが、来年こそはゴールしたい。

自宅

新居に引っ越すと共に固定ローラーを導入、冬期の練習不足を解消するのが狙い。ただ、1歳の息子がチェーンを鷲づかみにするので運用は難しい……。分譲マンションだが、自転車を各住戸まで持っていく前提でエレベーターや共用スペースが設計されている物件なので気に入っている。

自転車レースを観戦しながら練習可能な夢の環境が実現

自転車の良さ

ペダルを踏めばどこまでも行ける、というシンプルな楽しさはもちろんあるが。自分の限界と否応無しに向き会わなければならない点も良い。無茶した所で巡航速度は上がらないし、速く走るためには練習するための時間、良いパーツを買うための収入が必要である。人生のリソース配分最適化を考える良い契機になった。

まとめ

みなさん自転車やりましょう。こちらからは以上です。

このエントリーをはてなブックマークに追加

2015-10-13

PyConJPで「アドネットワークのデータ解析チームを支える技術」という発表をしました

PyConJP 2015に登壇してきました。PyConJPについてはライブラリの使い方についてのセッションが多く、じゃあそれを使ってどうやって金を稼いでるか、というプロダクトベースの話は少なめな印象がありました。なので自分は逆張りで仕事の話をしてきました。
といっても発表で例に出した手法は "Online Advertising" で論文を探せば出てくるメジャーな奴なので、うちだけ特殊な事をやっている訳では無いです。 IMG_7502

発表資料


補遺

発表の途中で「BigQueryは謎テクノロジ」と端折ってしまいましたが、O'reillyのBigQuery本に数テラバイトを数秒でスキャンする方法は書かれています。興味がある方はそちらを参照してもらえれば。

Google BigQueryGoogle BigQuery
Jordan Tigani,Siddartha Naidu,Sky株式会社 玉川 竜司
オライリージャパン
売り上げランキング : 260401
Amazonで詳しく見る by AZlink

反省

カンファレンス1日目と2日目も発表資料の仕上げに追われていて、他の人の発表をしっかりと聴講する余裕が無かったのが本当に悔やまれる所。夜のパーティや二日目のコーヒーブレイクでも他の参加者と話す事もできず、ポスターセッションも行けずじまい。次回はカンファレンスが始まるまでに仕上げておくと心に固く誓ったのでした。

 Link



このエントリーをはてなブックマークに追加

2015-07-05

Thompson Samplingの実験的な評価

Thompson Samplingの実戦投入の参考になりそうなので次の論文を精読した。
内容はアルゴリズムの紹介、パラメータチューニングとそれによる損失の変化の実験結果、オンライン広告とリコメンドエンジンへの適用例など。以降、Thompson SamplingはTSと表記。

UCBとの比較

実験ではどのパターンも、途中まではUCBと同じ損失だが試行回数がある程度増えた所でTSの損失がUCBを下回っている。腕の数が多い程、報酬の差が小さい程UCBとパフォーマンスが分かれる点は後ろにずれる。
理由は書いてないが、UCBは一旦報酬が低い事がわかった腕も定期的に引くのでこういう結果になるのだろう。TSは報酬の期待値が収束した後は良い物しか引かない。

事後確率の調節

ベータ分布のパラメータaとbについてそれぞれα∈(2, 1, 0.5, 0.25) で割った値にした実験。0.25にすると傾向として損失は下がるが、損失が増える駄目なケースも増える。腕の性能評価を早めるので、誤って評価してしまうのだろう。バリアンスが高くなるので玄人向けのチューニングに見える。

報酬のフィードバック遅延に対するロバスト性

腕を引いて即時に報酬が観測できない事はシステムの制約上良くあるよね、という前段に続いて報酬のフィードバックが遅れた場合の実験結果が載っている。1000ステップ遅れた場合にUCBの損失は10倍近くなるのに対して、TSは6倍程度。TSの方が性能良いという結果。

オンライン広告におけるCTR予測との組みあわせ

Contextual Banditな話。CTRの事後確率をベータ分布では無く、ロジスティック回帰で平均を求め、ラプラス近似でガウス分布にした物を利用する。PRMLでラプラス近似を見た時は何に使うのか理解できなかったが、なるほど納得。

この論文に書いてない事

腕の報酬が変化するパターンの実験はしていない。広告配信では腕の性能が変化する事は考慮しなければならないので気になる所。報酬の変化を捉えるにはUCBの様に一旦低く評価された腕を再度引くしかけが必要になるので、期待報酬の分散を増やすか利用する観測値は直近の物に限定する、といった工夫が必要だろう。

このエントリーをはてなブックマークに追加

2015-06-26

ISAC TokyoでGPS/GNSS高精度測位の世界を垣間みた

今年のInternational Space Apps Challenge、グローバル審査の結果まで出たので顛末を書きます。結果は東京ローカル審査で優秀賞を頂いて、グローバル審査に進んだものの敗退といった所です。

International Space Apps Challenge(ISAC)とは

今年で4回目となるNASA主催の国際ハッカソンで、今年は世界133ヶ所で行なわれた模様。地球観測衛星から得られたデータを有効活用したり、世界規模の問題解決に取り組むといったテーマが与えられた中でチームを作りハードウェアやWebサービスの開発に取り組みます。

何を作ったか

SENSOR YOURSELF テーマを選び、屋外用の安価な高精度測位モジュールを作りました。モチベーションとしては既存の機器、例えば農機や建設機器のオプションとして売られている物があまりに高いため、高精度測位の普及が進んでいない。安価に実現可能な事と利用法を示せれば、廉価品を作るメーカーが出てきたりスマホだけで利用可能になって嬉しいのでは、という話です。農機になぜ高精度測位が必要かと言うと、例えば植えた稲を踏まないように自動運転させるためだそうで。

作りとしてはRaspberry Pi + GPS/GNSS受信機 + ケースです。ネットワーク必須なので、そこはテザリングで。スマホ内蔵の受信機の信号が直接使えるようになればラズパイも受信機も不要なのですが、iOSはCore Location経由でしか位置情報にアクセスできないので不可能なのと、Androidでも特定の機種に限られるのでこの形に。測位モジュールを作るだけでなく、有効なユースケースを提示する予定だったのですがハッカソン中に良い案が思いつかず。

これでどの程度の精度が出るかと言うと、単独測位の誤差が10メートル程度なのに対して誤差1 ~ 10cm程度です。デモページで取得したデータをプロットしているのですが、1cm程度に収まっています。ただ、これは建物の屋上の受信機で測定したデータなので(Height 117m)地上でやるよりも精度が出ています。
デモのスクリーンショット

READMEに載せた構成図
 +------------------+          +--------------------------------------+                            
 | Smartphone       |          | RaspberryPi                          |                            
 | +--------------+ | Location |  +------------+       +-----------+  |           +-----------------+
 | |   Browser    <----(HTTP)-----+ API Server <-(TCP)-+  RTKLIB   <----(Serial)--+ GNSS/GPS Device |
 | +------^-------+ |          |  +------------+       +-----^-----+  |           +-----------------+
 +--------|---------+          +-----------------------------|--------+                            
          |                                                  |                                     
          |                                               (NTRIP)                                     
          |                                                  |                                     
 +--------+---------+                               +--------+---------+                           
 | Web Server       |                               |  Base Station    |                           
 | HTML, JavaScript |                               |                  |                           
 +------------------+                               +------------------+                           
 Serve User Interface
Submitしたプロジェクトページ

最近はGPSでは無くGNSSと呼ぶらしい

今回のハッカソンで勉強になったのがGPS周りの技術。基本的な所ではGPSは米軍の衛星測位システムを指すが、ロシアのGLONASSやEUが稼動準備中のGalileoと利用可能な衛星測位システムは複数存在するため、今はGPSでは無く総称のGNSS(Global Navigation Satellite System)を使うのが時流らしい。といっても文献ではGPS/GNSSと併記しているパターンが一番多かったので、この記事でも併記にしました。

単独測位と高精度測位

スマートフォンでWi-Fiを切った状態での位置情報の正確さ、地図上で小さな通りだと一本ずれてしまう事は良くあります。受信機単体測位による精度は良くて10メートルとされているので、普段の感覚通りですね。この時、4つ以上の衛星との距離を元に3次元空間上の位置を計算しています。距離は信号伝播時間と光速度の積で求まります。3点からの距離さえわかれば位置は計算できるはずですが、受信機の時計補正のために4台目の衛星を必要としています。

この精度を上げる方法がいくつかあり、DGPSだとメートル単位、RTKだとセンチメートル単位での測位が可能になります。
  • DGPS (Differential GPS) 
  • PPP (Precise Point Positioning) 
  • RTK (Real Time Kinematic GPS)
今回のハッカソンでは海洋大のRTK基準局が使えるとの事で、RTKにしました。ただ、RTKの基準局の有効範囲は15km程度なので、基準局がどこにでもある時代にならないとこの方式は普及しないのでは、とも思ったり。
出展: RTK-GPSの原理と応用

渋谷とIngressとマルチパス

これらの高精度測位が使えれば、常に数十メートルのドリフト(位置のずれ)によって苦戦を強いられている渋谷IngressエージェントもUltra Strikeが有効活用できそうですが、ビルの谷間では上手くいきません。前述の高精度測位手法は、誤差を発生させる要因の1つである反射波(マルチパス)をケアしないからです。ハッカソン会場は秋葉原だったのですが、地上においてRTK方式では一度も位置がFixしませんでした。RTK方式が利用できるのは、空が開けた場所に限定される様です。

ハッカソン参加で得られた物

1年目は二行軌道要素(TLE)について知り、2年目は太陽光パネルの発電量計算モデルと電力買い取りサービス、Modis(地球観測システム)のデータとGIS。今年はGPS/GNSS。普段の仕事では組まない人とプロジェクトをやると新しい知見が得られて面白い。あと、気になってるけど触った事のないライブラリを試すのに良くて、今年はRxJSを投入。毎年何らかの学びがあるので楽しみにしています。

状況です。

関連プロジェクト

どちらも専用ハードで小型化を実現していて凄い。海外だとRTK基準局がそこら中にあるのか、もしくは基準局を自分で立てて運用できるGPSおじさん向けのプロダクトかもしれない。

参考文献

坂井 丈泰: GPS/GNSSの基礎知識
土井下健治, 村本英一, 神田俊彦: 建設機械への ICT 応用
浪江 宏宗: RTK-GPSの原理と応用
安田 明生 et al, 東京湾仮想RTK基準局ネットワーク

このエントリーをはてなブックマークに追加

2015-06-12

弱くてニューゲームしてアドテクエンジニアになりました

近況。2015年からアドネットワークのデータサイエンスチームにおります。前の部署ではメディア寄りの所でモバイルアプリの開発をしていたので、ほぼ転職に近い状態です。アドネットワークなにそれという方向けの説明としては、広告主と広告枠をまとめていい感じにディスプレイ広告を配信するシステムだと思ってもらえれば。

現在のミッション

データ分析や広告配信アルゴリズムの改良というアプローチでアドネットワークの収益改善に取り組むのがミッションです。会社ブログにMortal Multi-Armed Banditsの記事を書いた頃は多腕バンディットアルゴリズムの調査や実装をしていました。

データサイエンスといっても、いきなり機械学習を使った仕組みをプロダクションに投入できるかというと全くそんな事は無く、ログの収集と解析基盤を構築する所からでした。まっさらなAWSアカウントでCloudFormationテンプレートを書いて、VPCやサブネットを切っていたのが2月頃。その間に並行して勉強も進めました。

弱くてニューゲーム

データ解析系の業務に興味はあったものの、経験は無く実力としては新卒以下。2011年に画像認識の論文で機械学習を知り、興味本位で勉強会に参加してみたもののPRML(後述)はニューラルネットが理解できず挫折し、Andrew Ng先生のオンラインコースは課題提出が間にあわず途中で落第する始末。それから3年経って仕事になるとは、どう転ぶかわからない物です。

解析系スキルとドメイン知識は1から修得しないといけない一方、新しいデータを得るために広告表示のJavaScriptコードに機能を追加したり、Scalaで動いてる配信サーバーにも手を入れたりします。このあたりは人に頼まなくても自分でやれてしまうので楽。広告表示のJavaScriptは流行りのライブラリなど一切使わないので、ECMAScript標準とDOMをおさえておけば書ける、元DOM職人としては血が騒ぎますね。

面白い所

扱うデータの量が今まで(そこまでヒットしていないサービス開発)とは桁が違う。マルチコアを使い切る、かつサーバーを横に並べればスケールするプログラムを書き、c4.8xlarge複数台並べてログの処理をするのは爽快。しかしそれも最初だけなので、新鮮な気持は忘れないでおきたいですね。

あとは論文が読める点。アドテク業界はプレイヤー同士がeCPMを高めるために金と計算機とアルゴリズムで殴り合いをしている様な物ですが、論文という形で成果が発表されるので、読みまくれば勉強になる。読むだけなのも何なので数年以内には書いてみたい。

読ん(だ|でいる)本をいくつか

ザ・アドテクノロジー~データマーケティングの基礎からアトリビューションの概念までザ・アドテクノロジー~データマーケティングの基礎からアトリビューションの概念まで
菅原健一,有園雄一,岡田吉弘,杉原剛
翔泳社
Amazonで詳しく見る by AZlink
アドテクノロジー プロフェッショナル養成読本~デジタルマーケティング時代の広告効果を最適化! Software Design Plusアドテクノロジー プロフェッショナル養成読本~デジタルマーケティング時代の広告効果を最適化! Software Design Plus
養成読本編集部
技術評論社
Amazonで詳しく見る by AZlink
どちらもネット広告の進化の歴史と要素技術をざっと見わたせる。THE AD TECHNOLOGYは広告全体におけるオンライン広告の位置づけ、マス広告・オフライン広告との関連にページが割かれているので広告業界初心者には良い。

プロフェッショナル養成読本シリーズは、どちらかと言えばマーケッター向けの内容だった。表紙がスーツマンなのはそういう事か。クロスデバイスターゲティング等のスマホ回りの仕掛けは載っていないので、スマートフォン向けの広告をやっている人は物足りないかも。

パターン認識と機械学習 上パターン認識と機械学習 上
C.M. ビショップ,元田 浩,栗田 多喜夫,樋口 知之,松本 裕治,村田 昇
丸善出版
Amazonで詳しく見る by AZlink
基礎を身につけるための定番の1冊。3年前は挫折したものの、今回はそうも言ってられない。とにかく数式をコードに落としてしまえば、動作もするし容易に理解できるのでなるべく手を動かす事にしている。最近だとQiitaに投稿したベイズ線形回帰によるパラメータ分布の収束がそれ。毎週の社内勉強会で読み進めているので2年もすれば下巻まで終りそう。これのおかげか、論文に書いてある事が理解できない事案が減りました。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
久保 拓弥
岩波書店
Amazonで詳しく見る by AZlink
PRMLはガウス分布が多めですが、こちらは何でもかんでもガウスで当てはめが上手く行く訳がない、という主張。ポワソン分布, 二項分布, ガンマ分布, etc. といった確率分布と適用例を紹介し、いかにデータの特徴をとらえた確率分布・モデルを選ぶかという事に主眼がおかれている。機械的に素性選択をする手法ばかり見た後に読んだので、自分にとってはカウンターカルチャー的な存在。

深層学習 (機械学習プロフェッショナルシリーズ)深層学習 (機械学習プロフェッショナルシリーズ)
岡谷 貴之
Amazonで詳しく見る by AZlink
次の輪読会でいきなり発表担当。既に背水の陣。

UDEMY ベイズ推定とグラフィカルモデル:コンピュータビジョン基礎1
こちらは進捗51%。Computer Vision向けの内容だがとにかくわかりやすい。PRMLで詰まったら立ち戻るのに良い。

まとめ

レベル0から出直しといいつつも、稼いでなんぼのエンジニア。事業にコミットしつつ、先生きのこるための鍛錬に励みます。

このエントリーをはてなブックマークに追加

2015-05-16

データ分析コンペについてLTをした

Powerd by enjo-generator
といった内容でした。

まだまだ修行あるのみ。週末はkddcup2015に向けて助走を開始します。


このエントリーをはてなブックマークに追加

2015-05-10

Google Cloud StorageとBigQueryで困ったら、まずはStack Overflowに書くと良さげ

GCPでサポートパッケージを購入していない人向けの話。

困ったらどうするか

GCPのサポートページを見ると、課金しないで利用できるのはドキュメントとコミュニティフォーラムとある。コミュニティフォーラムとしてGoogleグループとStack Overflowへのリンクが載っている、がGoogleグループはStack Overflowへ移行したとあり、現状機能はしていない。さらにBigQueryについてはGoogle Codeにissue trackerがある。

まとめると窓口は

Stack Overflowが一番反応が早い

自分が投稿した質問
こういう質問を投稿すると数日で「まだ現象が再現するならメール頂戴」みたいなレスが付く。他の質問を見ても、解答もしくはなんらなかのコメントが付いている。

Stack Overflowよりも先にGoolge Codeの方にも投稿したがこちらは3週間経っても反応が無い。よって、まずはStack Overflowに投稿するのが良いだろう。Stack Overflowは日本語版も出ているが、GCPについてのやりとりは活発では無い様子。

まとめ

Stack Overflowで解決しなかったら金での解決(サポートパッケージ購入)を検討しよう。

このエントリーをはてなブックマークに追加

2015-05-09

情報セキュリティスペシャリスト試験に合格していた

いつの間にか合格通知が届いていた、昨年10月の試験にパスしていた模様。試験自体については午後試験が簡単すぎるという印象のみが強く残っているが、思い出しながら書く。

受験動機

テキストがハイパフォーマンス ブラウザネットワーキングの復習になりそうだったので。ついでに試験も受けた。

テキスト

テキストは情報処理教科書の2015年版を選び、試験前に2週間かけて読んだ。試験に出る問題に限定したくなかったので一番ボリュームのある奴で正解。TCP、UDPの攻撃パターンやTLS、PKI・証明書チェーンについて一通りページがあったのでそこは満足。ISMS認証取得の実務について事細かく載っているが、生憎この内容が役に立つ予定は無い。

メールについては最近使っていなかったので補強になった。ここ3年はモバイルアプリケーションの開発をしていたため、メールを送信する代わりにPush通知とSMSばかり使っていた。

出題問題について

XSS, SQLインジェクションといった頻出問題は普段Webアプリケーションのコードを書いていれば当然解ける。しかし問題として登場するプログラムが酷すぎて読むだけで激しいストレスを受ける。例えばJava Servlet中にインラインでHTMLが書かれており、せめてJSPに分離しろという気分に。この様なコードが存在する現場はセキュリティ以前の問題が山積みだろう。関係ないけど出題されるJavaScriptコードでAngularJSが使われていたら2014年ぽくて良いのにな、などと考えていた。

記述問題のある午後試験、マネジメントについては問題文を良く読めば解ける問題、セキュアプログラミングについては常識レベルに留まっていたのでもっと難しくて良いと感じた。Webブラウザがhttps通信で安全性を確保している仕組みを5000字程度で述べよ、みたいな論述問題があっても面白いのではないか。

まとめ

テキストは興味深く読めて良かった。次は統計検定とLPICが気になるが、LPICは受験料が高いのでテキストで自習のみになりそう。

このエントリーをはてなブックマークに追加

2015-03-04

「データ解析のための統計モデリング入門」をPythonでトレースする (3章)

「データ解析のための統計モデリング入門」、もっと早くに取り組んでおけば良かった。データ分析業のアプローチに必要な観点を学べる良テキストだ。しかし本文中はRが使われているのでPythonで一通り書き直して読み進めてみる。
GLMフィッティングの際の勾配降下、確率的勾配降下法とかいろいろあったと記憶しているが、最もナイーブな実装しか書けなかった。修行が足りない。しかしPythonのstatsmodelsは今回初めて使ったが便利すぎる。


3.3 統計モデリングの前にデータを図示する

In [5]:
df = pd.read_csv('./data3a.csv')
In [12]:
def draw_scatter(df):
    plt.scatter(df[df.f=='C'].x, df[df.f=='C'].y, color='r')
    plt.scatter(df[df.f=='T'].x, df[df.f=='T'].y, color='b')
    plt.legend([u'f == C', u'f == T(施肥)'], loc='upper left')
In [13]:
draw_scatter(df)
plt.title(u'図3.2')
Out[13]:
<matplotlib.text.Text at 0x10d1812d0>
In [15]:
df.groupby(df.f).boxplot(column='y')
plt.ylim(0, 15)
Out[15]:
(0, 15)

3.4.1 線形予測と対数リンク関数

In [48]:
X = np.linspace(-4, 5, 100)
plt.plot(X, np.exp(-2 - 0.8*X), 'k--')
plt.plot(X, np.exp(-1 + 0.4*X), 'k-')
plt.legend(['$\{beta_1, beta_2\} = \{-2, -0.8\}$', '$\{beta_1, beta_2\} = \{-1, 0.4\}$'])
plt.title(u'図3.4')
Out[48]:
<matplotlib.text.Text at 0x113c84450>

3.4.2 あてはめと、あてはまりの良さ

In [49]:
def log_likelifood(b1, b2, x, y):
    # ポワソン回帰の対数尤度を返す
    ret = 0
    for i in range(x.size):
        ret += y[i] * (b1 + b2*x[i]) - np.exp(b1 + b2*x[i]) - sum(np.log(range(1, y[i])))
    return ret
独自で対数尤度最大化を実装してみる
In [59]:
def climb_b1(b1, b2, x, y):
    # b2固定でb1を動かした時の最大の対数尤度を返す
    sign = +1
    alpha = 0.01

    logL = log_likelifood(b1, b2, x, y)
    
    while True:
        b1 += alpha * sign
        logL_n = log_likelifood(b1, b2, x, y)
        if np.abs(logL - logL_n) < 0.00001:
            break
        if logL_n < logL:
            sign = sign * -1
            alpha = alpha * 0.5
        else:
            alpha = alpha * 1.2
        logL = logL_n
    return logL, b1

def hill_climb(alpha, x, y):
    b1 = 1
    b2 = 0.1
    sign = +1
    
    max_logL, b1 = climb_b1(b1, b2, x, y)
    
    while True:
        b2 += alpha * sign
        # b2を動かした時の対数尤度を出す
        max_logL_n, b1 = climb_b1(b1, b2, x, y)
        # 変化しなくなったら終了
        if np.abs(max_logL_n - max_logL) < 0.00001:
            break
        if max_logL_n < max_logL:
            # 対数尤度が減ったら符号逆転
            sign = sign * -1
            alpha = alpha * 0.5
        else:
            alpha = alpha * 1.2            
        max_logL = max_logL_n
    print('Likelihood: %s, b1: %s, b2: %s' % (max_logL, b1, b2))
    return b1, b2

hill_climb(0.1, df.x.values, df.y.values)
Likelihood: -35.8338102468, b1: 1.29054140804, b2: 0.0757920975997
Out[59]:
(1.2905414080350648, 0.07579209759971002)
statsmodelsを使うパターン
In [60]:
import statsmodels.api as sm
from statsmodels.formula.api import glm

def calc_x_related_model(df):
    return glm('y ~ x', data=df, family=sm.families.Poisson(sm.families.links.log)).fit()
In [61]:
x_related_model = calc_x_related_model(df)
x_related_model.summary()
Out[61]:
Generalized Linear Model Regression Results
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -235.39
Date: Wed, 04 Mar 2015 Deviance: 84.993
Time: 23:17:23 Pearson chi2: 83.8
No. Iterations: 7
coef std err z P>|z| [95.0% Conf. Int.]
Intercept 1.2917 0.364 3.552 0.000 0.579 2.005
x 0.0757 0.036 2.125 0.034 0.006 0.145
In [62]:
import scipy.stats
def plot_fig36(model):
    X = np.linspace(-0.1, 1.6, 200)
    plt.plot(X, scipy.stats.norm.pdf(X, loc=model.params[1], scale=model.bse[1]))
    plt.plot(X, scipy.stats.norm.pdf(X, loc=model.params.Intercept, scale=model.bse.Intercept))
    plt.title(u'図3.6')
    
plot_fig36(x_related_model)
In [63]:
def plot_likelifood():
    X = np.linspace(0.05, 0.10, 1000)
    Y = []
    for x in X:
        Y.append(log_likelifood(1.2926, x, df.x.values, df.y.values))
    plt.plot(X, Y)
plot_likelifood()
plt.title(u'尤度が最大になる付近のグラフは正規分布に近い?')
plt.xlabel('b1')
plt.ylabel(u'対数尤度')
Out[63]:
<matplotlib.text.Text at 0x114eaf110>
対数尤度の評価
In [64]:
x_related_model.llf
Out[64]:
-235.38625076986077
AICで評価
In [65]:
x_related_model.aic
Out[65]:
474.77250153972153

3.4.3 ポアソン回帰モデルによる予測

In [66]:
def draw_estimate_x(model):
    X = np.linspace(7, 13, 100)
    Y = model.predict({'x': X})
    plt.plot(X, Y, 'k-')
In [67]:
draw_scatter(df)
draw_estimate_x(x_related_model)
plt.xlim(7, 13)
plt.title(u'体サイズ$xi$を使った統計モデルによる予測')
Out[67]:
<matplotlib.text.Text at 0x1175f39d0>

3.5 説明変数が因子型の統計モデル

In [68]:
def calc_f_related_model(df):
    return glm('y ~ f', data=df, family=sm.families.Poisson(sm.families.links.log)).fit()
In [69]:
f_related_model = calc_f_related_model(df)
f_related_model.summary()
Out[69]:
Generalized Linear Model Regression Results
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -237.63
Date: Wed, 04 Mar 2015 Deviance: 89.475
Time: 23:17:31 Pearson chi2: 87.1
No. Iterations: 7
coef std err z P>|z| [95.0% Conf. Int.]
Intercept 2.0516 0.051 40.463 0.000 1.952 2.151
f[T.T] 0.0128 0.071 0.179 0.858 -0.127 0.153
In [70]:
def draw_estimate_f(model):
    X = np.linspace(7, 14, 1000)
    Y = model.predict({'x': X, 'f': ['C']*len(X)})
    plt.plot(X, Y, 'r-')
    Y = model.predict({'x': X, 'f': ['T']*len(X)})
    plt.plot(X, Y, 'b-')
In [71]:
draw_scatter(df)
draw_estimate_f(f_related_model)
plt.xlim(7, 14)
plt.ylim(4, 10)
plt.title(u'施肥効果$fi$を使った統計モデルによる推定')
Out[71]:
<matplotlib.text.Text at 0x118a22710>

3.6 説明変数が数量型 + 因子型の統計モデル

In [72]:
def calc_x_f_related_model(df):
    return glm('y ~ x + f', data=df, family=sm.families.Poisson(sm.families.links.log)).fit()
In [73]:
x_f_related_model = calc_x_f_related_model(df)
x_f_related_model.summary()
Out[73]:
Generalized Linear Model Regression Results
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 97
Model Family: Poisson Df Model: 2
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -235.29
Date: Wed, 04 Mar 2015 Deviance: 84.808
Time: 23:17:42 Pearson chi2: 83.8
No. Iterations: 7
coef std err z P>|z| [95.0% Conf. Int.]
Intercept 1.2631 0.370 3.417 0.001 0.539 1.988
f[T.T] -0.0320 0.074 -0.430 0.667 -0.178 0.114
x 0.0801 0.037 2.162 0.031 0.007 0.153
In [74]:
def draw_estimate_x_f(model):
    X = np.linspace(7, 13, 100)
    Y = model.predict({'x': X, 'f': ['C']*len(X)})
    plt.plot(X, Y, 'r-')
    Y = model.predict({'x': X, 'f': ['T']*len(X)})
    plt.plot(X, Y, 'b-')
In [75]:
draw_scatter(df)
draw_estimate_x_f(x_f_related_model)
plt.xlim(7, 13)
plt.title(u'体サイズ$xi$と施肥効果$fi$を組みこんだ統計モデルによる推定')
Out[75]:
<matplotlib.text.Text at 0x10fd15ed0>
In [44]:
glm('y ~ x + f', data=df, family=sm.families.Poisson(sm.families.links.identity)).fit().summary()
Out[44]:
Generalized Linear Model Regression Results
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 97
Model Family: Poisson Df Model: 2
Link Function: identity Scale: 1.0
Method: IRLS Log-Likelihood: -235.16
Date: Wed, 04 Mar 2015 Deviance: 84.538
Time: 20:13:50 Pearson chi2: 83.6
No. Iterations: 7
coef std err z P>|z| [95.0% Conf. Int.]
Intercept 1.2671 2.843 0.446 0.656 -4.306 6.840
f[T.T] -0.2048 0.582 -0.352 0.725 -1.346 0.936
x 0.6606 0.290 2.281 0.023 0.093 1.228

このエントリーをはてなブックマークに追加