富士山の初冠雪日(2) トレンド

前回作ったモジュール(read_FirstSnowFuji)で初冠雪日を読み、年初からの日にち（1月１日を１とする）に変換したうえで、線形回帰を当てはめてみる。
以前やったように、rpy2を用いて、Rにデータを渡してやる。

import numpy as np
import matplotlib.pyplot as plt
import scikits.timeseries as ts
#module
from read_FirstSnowFuji import *

#for R
import rpy2.robjects as robjects

ayear,amonth,aday= read_FirstSnowFuji()
days=np.array([],dtype="i4")
for  n in range(ayear.size):
    td=ts.Date("D",year=ayear[n],month=amonth[n],day=aday[n])
    days=np.append(days,td.day_of_year)


#for R
ryear=robjects.FloatVector(ayear)
rdays=robjects.FloatVector(days)
robjects.globalEnv["ryear"] = ryear
robjects.globalEnv["rdays"] = rdays
r=robjects.r

scripts="""
png('FujiSnow_trend.png')
plot(ryear,rdays,xlab="Year",ylab="Day of year") # scatter plot
fm<-lm(rdays ~ ryear) # linear model
abline(fm) # plot regression line                 
dev.off()             
"""
r(scripts)
summary=r('sm<-summary(fm)')
print summary

出力(print summary)は

Call:
lm(formula = rdays ~ ryear)

Residuals:
     Min       1Q   Median       3Q      Max
-54.2679  -5.8555   0.9491   8.1001  27.3063

Coefficients:
             Estimate Std. Error t value Pr(>|t|) 
(Intercept) 138.25153   71.89542   1.923   0.0570 .
ryear         0.06873    0.03683   1.866   0.0645 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.45 on 115 degrees of freedom
Multiple R-squared: 0.0294,    Adjusted R-squared: 0.02096
F-statistic: 3.484 on 1 and 115 DF,  p-value: 0.06453

係数が0.06873/yearだから、100年に一週間ほどのトレンド（初冠雪日が遅くなる）となる。ただし、p値が約0.0645であるから危険率5%でトレンドがない可能性を否定できない。データのばらつきが大きいことから、さもありなんといったところである。