2015年11月27日金曜日

sampling

あるとき、統計分析のため、
重みに応じたリサンプリングをする必要ができました。


floatのベクトルvがあったとすると

p=v/sum(v)

を確率ベクトルとして、復元抽出するコードを考えてみました。
探しては見たのですが、すぐには見つかりませんでした。
そこで、したのような処理で実現してみました。

import numpy as np

n = 100
def resamp( v, n ):
   vsum = np.sum(v)
   vcum = np.cumsum(v)
   indices = np.array([ np.argmax(1.0/(vcum - r)) for r in np.random.uniform(0,vsum,n) ])
   return indices


確率ベクトルvと出力サンプル数nを入力し、インデックスのnumpy.arrayを返す関数です。

備忘録のため、ブログに載せてみましたが、
関数一つでできる方法がないのだろうか。



2015年11月25日水曜日

regression OLS

pythonで回帰分析(regression)を行うことができます

回帰分析は古くから知られており、汎用的で世間で一番用いられている統計分析法かもしれません。
回帰分析というのは、簡単に言うならば、説明変数から目的変数をダイレクトに計算する線形式を決めることと言えます。
もちろん、その間に因果関係がなさそうなデータでも、入力(説明)変数と出力(目的)変数のデータさえそろえば、計算できます。 なので、風の風速と桶(おけ)店の振り上げの関係式も分析できます。

推定性能についてはR二乗値などの値で評価します。 ここでは「通常の方法」についてpythonで行う方法についてご紹介します。 OLS法というのはOrdinary Least Square(普通最小二乗法)です。 回帰分析は最小二乗法で「行うべきもの」と考えている人もいるかもしれませんが、 他にもいくつか選択肢はあります。 例えば、重み付き最小二乗とか、ロバスト推定とか。 なので、一番使われている「通常の方法」がOrdinary Least Squareです。


import numpy
import statsmodels.formula.api as sfa
d = {'x':numpy.random.randn(10,3), 'y':numpy.random.randn(10) }
o = sfa.ols( 'y~x', data=d )
r = o.fit()
print r.summary()


もし、上のコードでstatsmodelsが無いと言われたら、 easy_installで入れてください。 結果の表示は下のようになりました。
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.403
Model:                            OLS   Adj. R-squared:                  0.104
Method:                 Least Squares   F-statistic:                     1.348
Date:                Sat, 14 Nov 2015   Prob (F-statistic):              0.345
Time:                        23:01:08   Log-Likelihood:                -13.518
No. Observations:                  10   AIC:                             35.04
Df Residuals:                       6   BIC:                             36.25
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept      1.2127      0.745      1.628      0.155        -0.610     3.035
x[0]           0.6436      0.379      1.697      0.141        -0.284     1.572
x[1]          -1.1225      0.837     -1.341      0.229        -3.171     0.926
x[2]          -1.2621      1.742     -0.724      0.496        -5.526     3.001
==============================================================================
Omnibus:                        0.620   Durbin-Watson:                   1.691
Prob(Omnibus):                  0.733   Jarque-Bera (JB):                0.417
Skew:                          -0.425   Prob(JB):                        0.812
Kurtosis:                       2.472   Cond. No.                         8.32
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
 specified.
出力がたくさんありすぎて、わからなくなってしまうかもしれません。 真ん中の少し下辺りに係数(x[*] coef)や切片(Intercept coef)の推定値が表示されています。
あとはAIC・BICや基本統計量が表示されています。 AIC・BICは値が低いほうがよりデータに近いことを示す情報量規準というものです。 それぞれの統計量についてはまた回を改めて説明します。 回帰分析関数クラスのolsへ入力する'y~x'というのはGNU Rと同様の書式です。 つまり、この場合、yが目的変数で右側のxが説明変数です。

回帰分析や検定については他をググってみてください。 もっと詳しく丁寧に書いてあるページがたくさんあります。

2015年11月14日土曜日

python 科学計算関係のライブラリWindows版


pythonは多くのライブラリをimportすることが多いです。
一つのことをするにしてもnumpyやscipy、またグラフを描くならmatplotlibとか。

PIPがうまく動作しない場合も多く、Windowsの場合インストーラを使わなければならないばあいが多いような気がします。 Linuxなら該当するパッケージ名さえわかれば、aptとpipでたいてい何とかなります。 Linxuの場合の話はまた今度します。
それぞれ、pipやeasy_installでインストールできるものもあれば、できないものもありややこしいので、現状でのインストール方法をまとめてみました(2015年10月現在)。

numpyやscipyはeasy_installやpipを使うとlapackなどのライブラリと接続できないなどというエラーがでて面倒なので、おとなしくサイトからsuperpackというものをダウンロードしてインストールしました。

ウィンドウズでバイナリをダウンロードしてくる場合はpythonのバージョン、32/64ビットの違いに注意してください。
pandasはeasy_installでうまくいかず、pipでインストールできました。 なぜかeasy_installではvcvarsall.batが無いと怒られます。そのときのエラーをコピペしておきます。
 UPDATING build\lib.win32-2.7\pandas/_version.py
 set build\lib.win32-2.7\pandas/_version.py to '0.17.0'
 error: Setup script exited with error: Unable to find vcvarsall.bat

統計分析のstatsmodelsはeasy_installでインストールしました。
ただし、scipyやpandas、patsyを使っているようなので、それぞれをインストールしてからにしましょう。
patsyはpipでできました。

そのほかにもいくつか科学技術計算ライブラリがありますが、 長くなりそうなので、残りはまた今度ご紹介します。