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

前回の続きとして、私自身あまりよく理解しているとは言えないが、ココを参考にして、ロバスト推定でも線形回帰をとってみた。

参考：
ロバスト推定（統計数理研究所）
http://www.ism.ac.jp/~fujisawa/research/robust.html
ロバスト推定（朱鷺の社）
http://ibisforest.org/index.php?%E3%83%AD%E3%83%90%E3%82%B9%E3%83%88%E6%8E%A8%E5%AE%9A

Rのパッケージrobustbaseを使うために

install.packages("robustbase")

を行っておく。
以下がプログラム。

import numpy as np
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="""
library("robustbase")   #import robustbase
png('FujiSnow_trend_robust.png')
plot(ryear,rdays,xlab="Year",ylab="Day of year") #scatter plot 
fm<-lmrob(rdays ~ ryear)  # robust linear model
abline(fm)  #plot regression line
abline(lm(rdays ~ ryear),lty=2) # usual linear model
dev.off()
"""
r(scripts)

# summary
summary=r('sm<-summary(fm)')
print summary

以下が図。点線が前回の回帰直線、実線がrobustbaseで得られた直線である。

 

出力(print summary)は

Call:
lmrob(formula = rdays ~ ryear)

Weighted Residuals:
    Min      1Q  Median      3Q     Max
-57.989  -7.872  -0.125   5.986  25.291

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 76.10515   51.85081   1.468 0.144896   
ryear        0.10154    0.02673   3.798 0.000234 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Robust residual standard error: 10.03
Convergence in 12 IRWLS iterations

Robustness weights:
 2 observations c(21,115) are outliers with |weight| <= 0.00057 ( < 0.00085);
 13 weights are ~= 1. The remaining 102 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.00365 0.88600 0.94480 0.88410 0.98680 0.99880
Algorithmic parameters:
tuning.chi         bb tuning.psi refine.tol    rel.tol
 1.5476400  0.5000000  4.6850610  0.0000001  0.0000001
 nResample     max.it     groups    n.group   best.r.s   k.fast.s      k.max
       500         50          5        400          2          1        200
 trace.lev compute.rd
         0          0
seed : int(0) 

これによると回帰係数は0.10154（100年で約10日冠雪日が遅くなる）で、しかもp値が非常に小さく統計的に有意である。
ただし、c(21,115)が外れ値、すなわち１９１４年（ayear[21-1]=1914)と２００８年（ayear[115-1]=2008）が外れ値になるのは良いとして、weightが１に近い点は13だけで、あとの外れ値を除く１０２点はweightが１以下に落とされている。

これを見るために、回帰直線の係数、weightなどの数値をpythonに戻し、プロットしてみる（同じことをRでもできるはずだが、pythonのmatplotlibのほうがなれているので）。
以下プログラム。
robust推定で得られた回帰直線を緑線で、
weightが1に近い点を青○で、
weightが0.5から0.9988の間にある点を黒点、
weightが0.00085から0.5の間にある点を赤点、
外れ値を赤xでプロットした。

#weights
weights=r('weights<-sm$weights')
weight=np.array(weights)

#plot
import matplotlib.pyplot as plt
plt.plot(ayear,ayear*a+b,'g-')
range=(weight>0.9988)
plt.plot(ayear[range],days[range],"bo")
range=np.logical_and(weight<=0.9988,weight>=0.5)
plt.plot(ayear[range],days[range],"k.")
range=np.logical_and(weight<0.5,weight>=0.00085)
plt.plot(ayear[range],days[range],"r.")
range=(weight<0.00085)
plt.plot(ayear[range],days[range],"rx")
plt.xlim(ayear[0],ayear[-1]) 
plt.xlabel("year")
plt.ylabel("day of year")
plt.savefig("FujiSnow_trend_robust_2.png")
plt.show()

ちなみに今年２０１０年の重みはweight[-1]=0.86999219949571671
このように多数の点の重みを落とすのは統計的に良くても、物理的な理由は自明ではない。
今回、この方法を用いたのは適当ではなかったかもしれない。

前回と今回の解析で、初冠雪日にトレンドがあるかは統計的にははっきりしなかったが、初冠雪日が遅くなっていっているかもしれないという作業仮説自体は興味深いものがある。これから、将来、どのような傾向があらわれてくるだろうか？

参考
異常気象を追う 初冠雪
http://www.bioweather.net/column/essay2/aw15.htm

縮小する富士山の永久凍土
http://www.fujisan-net.jp/data/article/1583.html