串本・浦神の潮位差と黒潮流路緯度 Scipyでregression

regressionはscipyでも行える。stats.linregress関数（実際にはmasked array用のstats.mstat.linregress)を使う。
参照　Cookbook/Linear Regression

#read module
from scipy import stats

from read_kuroshio_lat import read_kuroshio_lat
from read_kushimoto_uragami import read_kushimoto_uragami

#read data
kuroshio_lat=read_kuroshio_lat()
kushimoto_uragami=read_kushimoto_uragami()   

#calculate regression
(a_s,b_s,r,prob,stderr,c)=stats.mstats.linregress(kushimoto_uragami,kuroshio_lat)
print "slope=",a_s
print "intersect=",b_s
print "correlation=",r
print "p value=",prob
print "std. of residuals",stderr
print "(Var. Latitude)/(Var. Sea level)=",c

最後のcはマニュアルには無いが、ソースコードを見てみたところ(ipythonではstats.mstats.linregress??で見ることができる。scipy version 0.7.0）、それぞれの変数の分散の比になっていた。結果は

slope= 0.0666657
intersect= 31.9559306074
correlation= 0.618468
p value= 4.38452479663e-62
std. of residuals 0.62129592706
(Var. Latitude)/(Var. Sea level)= 0.011619

regressionの線も引いてみよう。

 

プログラムはmatplotlibを使って以下のようになる。

#plot
import matplotlib.pyplot as plt
import numpy as np

x=np.arange(-7.,35.,1.) # x data from -7 to 35 step 1
y=a_s*x+b_s         # regression line

plt.figure()
plt.plot(kushimoto_uragami,kuroshio_lat,'.') # plot data
plt.plot(x,y,'r-') # plot regression line
plt.xlabel("Kushimoto-Uragami")
plt.ylabel("Kuroshio Latitude")
plt.xlim((-7.,35.))
plt.ylim((30.,33.5))
plt.savefig("kurosho_reg_scipy.png")
plt.show()

今回はつかわなかったが、以下の情報もある。
Cookbook/Fit Statistics 

Scikits.statsmodels

multiple regression