tetsunosukeのnotebook

tetsunosukeのメモです

[pandas] pandas で 回帰分析

python pandas

Rで回帰分析 - tetsunosukeのnotebook を pandasでやってみた

特にファイルを読み込む部分がRっぽく書ける

import pandas as pd
>>> data = pd.read_csv("2-1.csv")
>>> data.describe()
          degree       amount
count  12.000000    12.000000
mean   16.391667  1310.833333
std     7.607587   414.355323
min     5.100000   772.000000
25%     9.725000  1127.500000
50%    16.700000  1231.000000
75%    22.750000  1355.000000
max    27.500000  2389.000000

[8 rows x 2 columns]

あとはRのnlsのように... 

>>> model = pd.ols(y=data["amount"], x=data["degree"], intercept=True)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "c:\Python27\lib\site-packages\pandas\stats\interface.py", line 135, in o
ls
    return klass(**kwargs)
  File "c:\Python27\lib\site-packages\pandas\stats\ols.py", line 53, in __init__

    import scikits.statsmodels.api as sm
ImportError: No module named scikits.statsmodels.api


ols関数で回帰分析しようとしたら落ちた！

追ってみると、Patsy が入っていなかったのでインストールして再実行

>>> model

-------------------------Summary of Regression Analysis-------------------------


Formula: Y ~ <x> + <intercept>

Number of Observations:         12
Number of Degrees of Freedom:   2

R-squared:         0.6119
Adj R-squared:     0.5731

Rmse:            270.7228

F-stat (1, 10):    15.7685, p-value:     0.0026

Degrees of Freedom: model 1, resid 10

-----------------------Summary of Estimated Coefficients------------------------

      Variable       Coef    Std Err     t-stat    p-value    CI 2.5%   CI 97.5%

--------------------------------------------------------------------------------

             x    42.6066    10.7296       3.97     0.0026    21.5766    63.6365

     intercept   612.4407   192.4569       3.18     0.0098   235.2251   989.6562

---------------------------------End of Summary---------------------------------

x, interceptの値、42.6066と612.4407を得ることが出来ました。


>>> model.summary_as_matrix
                 x   intercept
beta     42.606567  612.440687
p-value   0.002639    0.009782
std err  10.729551  192.456913
t-stat    3.970955    3.182222

[4 rows x 2 columns]

beta値のところに求めたい値が出ています。

R-Style

Rっぽくモデルを書くことができるので、それを試してみます。

>>> import statsmodels.formua.api as sm
>>> res = sm.ols(formula="amount ~ degree", data=csv).fit()
>>> res.summary2()
<class 'statsmodels.iolib.summary2.Summary'>
"""
               Results: Ordinary least squares
============================================================
Model:                 OLS     AIC:                 170.2930
Dependent Variable:    amount  BIC:                 171.2628
No. Observations:      12      Log-Likelihood:      -83.146
Df Model:              1       F-statistic:         15.77
Df Residuals:          10      Prob (F-statistic):  0.00264
R-squared:             0.612   Scale:               73291.
Adj. R-squared:        0.573
------------------------------------------------------------
           Coef.   Std.Err.   t    P>|t|   [0.025    0.975]
------------------------------------------------------------
Intercept 612.4407 192.4569 3.1822 0.0098 183.6200 1041.2614
degree     42.6066  10.7296 3.9710 0.0026  18.6996   66.5135
------------------------------------------------------------
Omnibus:             4.269      Durbin-Watson:         2.667
Prob(Omnibus):       0.118      Jarque-Bera (JB):      1.424
Skew:                0.700      Prob(JB):              0.491
Kurtosis:            3.943      Condition No.:         44
============================================================

Rでカイ二乗検定 を pandas/scipy で

python pandas

Rでカイ二乗検定 をPythonで行います。

カイ二乗検定を行うだけであればpandasは必要なさそうです。

Rで下記のように行っていたコードを、

chisq.test(
  c(10, 10, 10, 25, 1, 4)
)


scipy.stats.chi2_contingencyに渡すために、本来同数であれば全部10回ずつになる、という仮定を並べておきます。

>>> import scipy.stats
>>> obs = np.array(10,10,10,10,10,10], [10,10,10,25,1,4)
>>> res = scipy.stats.chi2_contingency(obs)
>>> res[0:2]
(16.363636363636363, 0.005879038879660156)

p値が0.005ということになりました。

同様に

10, 10, 10,  8,  8, 14

というデータであればp値は0.95311546950018988となりました。

検定の結果としては両方共Rで実施した場合と等しくなっています。

Rでこれと同じ数字を出すためには

> m = matrix(c(10,10,10,10,10,10, 10,10,10,25,1,4),ncol=6, byrow=T)
> m
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   10   10   10   10   10   10
[2,]   10   10   10   25    1    4
> chisq.test(m)

        Pearson's Chi-squared test

data:  m
X-squared = 16.3636, df = 5, p-value = 0.005879

このようにしないといけなかったのですね。

うどん、そば、ラーメンのケースでは

>>> scipy.stats.chi2_contingency(np.array([[1,7,3],[13,6,10]]))[0:2]
(7.5634710117468744, 0.022783116860954487)

となりRと結果が異なっています。これはR側で計算値が不正確、というメッセージのとおりだからかもしれません。

2013年を振り返る

会社のチームで、そもそも今年何したんだっけ、ってのを並べたらウンザリするほどろくな成果を出してない感じがあって、なんでだったんだろーねー。という感じ。

さて、全体としては

テスト。テスト。テスト。

結局チームでの在籍年数が一番長くなってしまったのですが、一向に技術的負債が減らないのをなんとか1日二時間残業してでもやる、というつもりで続けた。仕事がどうにもならなく忙しい時は出来なかったが、目標としてしっかり負い続けることができたのは良かった。

とはいえもとのコードに手を入れないと書けない部分などはそう簡単には手を付けられず、まあ、こうしてればよかったのに、ということのイメージが少しついたことぐらいがプラスか。

担当領域の変更

会社として新しいことをやっていくということで、リーダーからサブリーダーへ。サブリーダーというのはある意味性に合っていて、サポートをしていくということには変わりがなかったが、意志いれという部分で構想を持たなくなったのは自分としてよかったのかどうかは悩むところ。

また、新メディア、STANDS! で、マーケティング分野へ関わることになり、何かをはじめるにあたりデータの観点からどうサービスを構築するか？ということに強く関心を持ち、実業務においてもいろいろな経験することができた。これに関してはまだまだだけどとにかく経験だと思う。

データに関するお仕事

とにかくデータに関する本をいろいろ読んだ。データサイエンティストという名前のついた仕事にいろいろある中、自分の働き方にあったものがいくつか見えたので、それにデータサイエンティストという名前が不適切であったりそもそもそんな名前がクソだ的な話はいろいろあるが、職種であるかどうかも含めてここはやっていこうという気になった。

今年の冬休みの宿題はリーンアナリティクスを読むこと。
来年も事例にはいろいろ当たりたい。

健康について考え方を変える

しばらく続けていたトレーニングはやめて、もうすこし体に負担のかからないものにしようかと。
どうしても義務感から体調崩したりもしたので、無理なく続けられる何かに目を向けていきたい。
まあここは
夫婦でできるもの、というところへもシフトしたいということはあるのだけど。

趣味

今年は野球はたくさんの試合に出れるようにした。またまだ成績は良くならないけど、バッティングについてはいろいろ教えてもらってシーズン.260は打てたので及第点。守備が不安なのでなにかいい対策はないかなー。

リアル脱出ゲームについては少し慣れてきた。まあ慣れながらやっとクリアには二歩前、くらいのレベルだけど。

音ゲーは好きな曲しかしてなくて回数も少なかったからちっとも成長がないな。

2014年もやり過ぎない程度に続けたい。

2014年は

2013年少し感じたことだけど、自分が大事だと思ったことをなんとしてもやろう、という年にしたい。周りを見てもそういう人が上手くやっていたように見えた。

周りのことも大事だけど、そうでないことをしっかりやれるようにしたい。

データ運動会に行ってきた

データ運動会って？

こちらのイベントです。
http://atnd.org/event/E0019571

スポーツの秋ということで、ワタクシが見るのもやるのも大好きな野球のデータについてのアイデアソン・ハッカソンをやるということで参加してきました。

データ提供元はもちろんデータスタジアム様。

現状、プロ野球ファンは１０００万人、１６％くらい減っている状況で、テレビ放送が少なくなったりしている今、野球離れの進んだ若年層をターゲットにした野球を面白いと思わせ、野球市場を活性化するアプリを作ろうぜ、というテーマ。

データの可視化

データアーティスト様、アイオイクス様からも、データの可視化についての事例の紹介や、インフォグラフィックの注目度などのお話を頂きました。

データの見せ方3箇条
  • デザイン＜データ の関係は崩さない　：　何のデータかわかるようにする
  • データを入れすぎない
  • データは必ず信用できるものを使う

など、参考になりました。

アイデアソン

ラーニングプロセス　矢吹様、イノベーションファシリテータ・ワークショップデザイナーであるとのことで、マンダラートからはじめとして、他花受粉などのプロセスでとにかくアイデアを出しては交換し、最終的に投票で良いと思ったものをハッカソン用の元ネタとして使うということで。

個人的には、野球をどう身近にするか、ということで、イケメン度や好きになれそうな選手をレコメンドするアプリや、カードバトルなどで隠れたうんちく要素が発動（清原が泣き出して守れないとか、誰それは中学の先輩だから頭があがらないとか）するゲームとか、記録達成の日を予測してその日にスタジアムにいると記念品がもらえるとか、をアイデアとして出しました。

ハッカソン

自分が参加することになったチームは当初「監督よりも優れた采配ができるか？」というのをクイズ・投票形式にする、というアイデアでした。勝利に対しての寄与率みたいな数字がモデル化されているので、それが上下するかなどで算出できるかな？と思っていましたが、そもそも采配を当てたところで野球を見る人って楽しいのか？ということになり、再度議論した結果、一球速報に、それぞれの対戦でどちらが勝ちそうなのか？というのを判定するアプリにしました。勝ちそうな要因にパワプロなどである、特殊能力（初球・ランナー有・連発など）を可視化し、解説者がなぜ「今はバッターが有利」などというコメントをするのか、理解可能にする、という案に落ち着きました。

私以外のメンバーが統計系の大学院生だったので、モデルの作成は彼らに任せ、私はPHP+MySQL+JS+TwitterBootstrap での画面側を作りました。

下記のようなアプリができました。

f:id:kidd-number5:20131008121434p:plain

データスタジアム賞 を頂きました。実際にサービスを行っている一球速報をさらに改善できるアイデアであったということで、アイデアを練りなおして、ユーザ目線に何度も立ち返ったことが一つのポイントであったと思います。

その他のチームでは
  • ストライクゾーンの判定に偏りがあるかを調べ、審判が贔屓しているか？を調べるアプリ
  • 良妻判定　結婚した後の成績があがったかどうか、を、ピンチ・チャンス時の成績をもとにして判別
  • カープの応援であるスクワットをスマホセンサーで検知して、全国ランキングを作ったりするアプリ
  • 捕球した打球の重心からの距離、標準偏差などをもとに守備範囲を見える化する
  • 連打を浴びた、などの状態からパワプロのピヨり状態の推定
  • ストライクをボウリングと勘違いする人がいるのでそれなら対戦成績から打球の飛ぶゾーンを算出してボウリング
  • 選手の強さを定義してそれを倒した成績等から判定する隠れたお得意様を探すアプリ
  • ピタゴラス勝率などの数字から、ビール指数というビールがおいしく飲めそうな試合、タイミングを探るアプリ

上記概要だけだと工夫や面白さがちょっと伝わらないのですが、それぞれ独自の指標を作ったり、すでにある指標のビジュアル化、アプリとしてもスマホのセンサーを使ったりと、なるほどこういう着眼点で見せるというのは、同じデータを使っていてもやり方がいろいろあるのだなと、盗むべきところが多いイベントでした。
　

感想

メンバーによってアプリ屋がいなかったので、統計モデルの作成までにとどまったところや、初めてJavaScriptを書いた！みたいな人たちもいて、何かしらバランスを考えられたらよかっただろうな・・と思いつつも、それはそれで面白いかも・・・なんて。

しかし優勝したチームはアイデアソンの初日の後に飲みに行ったり、夜中まで実際に開発をしていたそうで、学生さんの若さと優秀さにただただ唖然としたし、同時に羨ましくも思いました。（おっと、学生さん以外も夜中までやっていた方がそのチームのリーダーでした。すごすぎる！）


また参加してみたいし、弊社にはスポーツ関連の写真が大量にありますので、これらを使ったハッカソンなどもご提供できると思います！


当日のイベントがNHKさんのニュースになっていました。
http://www3.nhk.or.jp/news/html/20131007/k10015088001000.html

メモ： Sass をWindowsに入れようとしたらgemで怒られたよ

デザイナーのデスクトップ（Windows）に、Sassを入れようとしたんだけど、

C:\work>gem install sass
ERROR:  While executing gem ... (Errno::EACCES)
    Permission denied - C:/Program Files/ruby-1.9.3/bin/sass.bat

インストールに失敗！！

したので、下記の手順で対応。

Ruby をインストールするときにできる、
Windowsのスタートメニュー→すべてのプログラム→ Ruby 1.9.3-p194
→ Start Command Prompt With Ruby を右クリック
→ 「管理者として実行」 すると、黒い画面になるので

C:\Windows\system32>gem install sass
Fetching: sass-3.2.10.gem (100%)
Successfully installed sass-3.2.10
1 gem installed
Installing ri documentation for sass-3.2.10...
Installing RDoc documentation for sass-3.2.10...

成功しました！

与えられた変数が数字かどうか調べる　のに ctype_digitを使って死んだ話

PHP

プログラムの受け入れテストをしていて、どうしても通らない箇所があって、こんな感じに書いてあったのです

if (ctype_digit($num)) {
    // 通ってくれない処理
}

まあもっと言うとこんな

$list = array(
    '1' => 10,
    '2' => 20,
);

foreach($list as $key => $value) {
    if (ctype_digit($key)) {
        // 通ってくれない処理
    }
}


というわけで

<?php
$list = array(
    '1' => 10,
    '2' => 20,
);

foreach($list as $key => $value) {
    var_dump($key);
    var_dump(ctype_digit($key));
}

int(1)
bool(false)
int(2)
bool(false)

ファッ？！

ああ、これはPHPの配列のキーがいつのまにかintegerになっちゃう問題なのねー。

というわけで ctype_digitは

この関数を活用するには string を渡さなければなりません。

ってやつだよ。

で、この$listはもともとこんな風にベタでは書いてなくて、実際には'98765432' みたいな値が入る。なので、動作テスト時にこれを書いた人は通ったんだろう。

で、なんで通るの。

<?php
var_dump(ctype_digit(98765432));

bool(true)

なんとっ

たとえば integer を渡すと、期待する結果にならない可能性があります。 HTML フォームに入力された整数値は、integer ではなく string 型で返されます。 マニュアルの 型 についての節を参照ください。

ああこれなのね

var_dump(ctype_digit(47));
var_dump(ctype_digit(48));

bool(false)
bool(true)

うわー。

<?php
for( $i = 1; $i <= 100; $i++) {
    if (ctype_digit($i) === TRUE) {
        echo 'T';
    } else {
        echo 'F';
    }

    if( $i % 10 === 0 ) {
        echo PHP_EOL;
    }
}

FFFFFFFFFF
FFFFFFFFFF
FFFFFFFFFF
FFFFFFFFFF
FFFFFFFTTT
TTTTTTTFFF
FFFFFFFFFF
FFFFFFFFFF
FFFFFFFFFF
FFFFFFFFFF

見事に48～57までナイジャー・モーガン。

var_dump(chr(47));
var_dump(chr(48));
var_dump(chr(49));
var_dump(chr(57));
var_dump(chr(58));

string(1) "/"
string(1) "0"
string(1) "1"
string(1) "9"
string(1) ":"

あーなるほどねー。

  • 128 から 255 までの整数値を渡すと、ひとつの文字の ASCII 値とみなします (負の値には 256 を足して、拡張 ASCII の範囲に収まるようにします)。 それ以外の整数値は、10 進数を含む文字列とみなします。


あぶねーなーもう

Web検索系のAPIって結構死んでるんだ・・・（BingWeb検索APIを試してみたよ）

RESTとかで使えるAPIを探せ

特定のキーワードで検索して上位3件のデータを取りたい！

みたいなことを思ってWeb検索系のAPIを探してみました。

Yahoo!


そういえばYahoo!の検索が有料になったんだっけとか思ったら、もう終了する予定の模様。

http://developer.yahoo.co.jp/webapi/search/

□提供終了日：
・2013年8月14日
※新規利用の申請については既に終了させていただいております。
※現在ご利用いただいているみなさまは、更新いただくことなく提供終了日までご利用いただけます。

Google

Googleさんは？

https://developers.google.com/web-search/?hl=ja

Note: The Google Web Search API has been officially deprecated as of November 1, 2010. It will continue to work as per our deprecation policy, but the number of requests you may make per day will be limited. Therefore, we encourage you to move to the new Custom Search API.

なんかもうCustomSearchAPI使って世界みたいね。両方共（まあ中身同じだから？）

それではCustomSearchAPIは？

Important: The Google Custom Search API requires the use of an API key, which you can get from the Google APIs console. The API provides 100 search queries per day for free. If you need more, you may sign up for billing in the console.

無料だと1日100回までなのかー。詰んだ。

そもそも自サイト内検索とかが主な目的みたいだしなんか違う気がする。

Bing

最初思いつかなかったけどこんなエンジンもありましたね

http://datamarket.azure.com/dataset/bing/search

ふむ。Windows Azure Marketplaceとやらに登録して利用することができるらしい。こっちは月に5000リクエストまでOKのようだ。Googleよりちょっとだけ多いレベル。

すんごいわかりにくいんだけどページ右側に Bing API Quick Start & Code　というのがある。

これを見てみると、.NETとかPHPのサンプルコードが載っていた。

あと、これもまたわかりにくいのだけど、（つかこのテキストリンクなのかよ）

このデータセットの参照
DataMarket サービス エクスプローラーを起動して、データセットを対話形式で参照できます。

を参照すると、実際にどういうURLを叩くとどういう結果が返ってくるか試せる。（トランザクション数がカウントされて減らされるので注意）

例えば

Query: 検索
Markter:ja-JP

を入れると、

こんなかんじのURLになるそうだ。
https://api.datamarket.azure.com/Bing/Search/v1/Web?Query=%27%E6%A4%9C%E7%B4%A2%27&Market=%27ja-JP%27

だが、このURLを直接叩くとベーシック認証を要求される。

このURLはベーシック認証で叩くものらしい。そのキーはといえば、

https://datamarket.azure.com/account/keys アカウントキーというところで取得できる。アカウントキーのくせに、ベーシック認証のID/パスワード両方ともこれらしい。なんだそら。

踏まえてプログラム

サンプルはそのまま動かなかったのでcurlバージョン。

<?php
// 検索したい文字列
$query = "検索";
// アカウントキー
$accountKey = '<key>';
// とりあえずWeb検索に限定。結果はJSONにする。「$format」なので注意。
$ServiceRootURL = 'https://api.datamarket.azure.com/Bing/Search/v1/Web?$format=json';
// ここはとりあえず検索キーワードだけ
// それぞれのパラメータをシングルクォートでくくる仕様
$queryParameters = array(
    'Query' => "'" . $query . "'",
);
$url = $ServiceRootURL. "&" . http_build_query($queryParameters);
$result = curl_get($url, $accountKey);

// 先頭の1件を出力
var_dump($result[0]);

function curl_get($url, $key) {
    $ch = curl_init();
    curl_setopt($ch, CURLOPT_URL, $url);
    curl_setopt($ch, CURLOPT_RETURNTRANSFER, 1);
    // ID/PWともにアカウントキーなので
    curl_setopt($ch, CURLOPT_USERPWD, $key. ":" . $key);
    curl_setopt($ch, CURLOPT_HTTP_VERSION, CURL_HTTP_VERSION_1_1);
    $buf = curl_exec($ch);
    curl_close($ch);
    $json = json_decode($buf);
    return $json->d->results;
}
結果
object(stdClass)#3 (6) {
  ["__metadata"]=>
  object(stdClass)#4 (2) {
    ["uri"]=>
    string(91) "https://api.datamarket.azure.com/Data.ashx/Bing/Search/v1/Web?Query='検索'&$skip=0&$top=1"
    ["type"]=>
    string(9) "WebResult"
  }
  ["ID"]=>
  string(36) "b8d565a7-0757-4603-8981-bc139c88e4b8"
  ["Title"]=>
  string(33) "郵便番号検索 - 日本郵便"
  ["Description"]=>
  string(270) "郵便番号検索はこちらから。地図、住所、郵便番号から郵便番号を検索できます。 English サイトマップ よくあるご質問・お問い合わせ ここからサイト内検索です 検索したい文字列を入力してください"
  ["DisplayUrl"]=>
  string(29) "www.post.japanpost.jp/zipcode"
  ["Url"]=>
  string(37) "http://www.post.japanpost.jp/zipcode/"
}