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
#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()
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
from read_NINO import read_NINO
#for R
import rpy2.robjects as robjects
data=read_NINO()
x=data["NINO3 ANOM"]
year=str(x.year[0])
month=str(x.month[0])
#data to R
nino3=robjects.FloatVector(x)
robjects.globalEnv["nino3"] = nino3
#r script
r=robjects.r
rscripts="""
nino3ts=ts(nino3,start=c(%(year)s,%(month)s),frequency=12)
png('NINO3_R_wavelet.png')
#import library wmtsa
library("wmtsa")
## calculate the CWT of the sunspots series using
## a Mexican hat wavelet (gaussian2)
a<-wavCWT(nino3ts)
## print the result
print(a)
## plot an image of the modulus of the CWT and the
## time series
plot(a,series=TRUE)
dev.off()
"""
r(rscripts %locals())
Continuous Wavelet Transform of nino3ts
---------------------------------------
Wavelet : Mexican Hat (Gaussian, second derivative)
Wavelet variance : 1
Length of series : 724
Sampling interval : 0.08333333
Number of scales : 73
Range of scales : 0.0833333333333333 to 60.3333333333333
#coding: utf-8
from read_NINO import read_NINO
#for R
import rpy2.robjects as robjects
data=read_NINO()
x=data["NINO3 ANOM"]
year=str(x.year[0])
month=str(x.month[0])
#data to R
nino3=robjects.FloatVector(x)
robjects.globalEnv["nino3"] = nino3
#r script
r=robjects.r
rscripts="""
nino3ts=ts(nino3,start=c(%(year)s,%(month)s),frequency=12)
png('NINO3_R_acf.png')
op<- par(mfrow=c(1,3)) # 描画設定：1行3列画面表示
acf(nino3ts,lag.max=72,xlab="lag (year)") # 自己相関プロット
acf(nino3ts,type="covariance",lag.max=72,xlab="lag (year)") #自己共分散プロット
acf(nino3ts,type="partial",lag.max=72,xlab="lag (year)") #偏相関プロット
par(op) #作業前の描画設定に戻す
dev.off()
"""
r(rscripts %locals())
from read_NINO import read_NINO
#for R
import rpy2.robjects as robjects
data=read_NINO()
x=data["NINO3 ANOM"]
year=str(x.year[0])
month=str(x.month[0])
#data to R
nino3=robjects.FloatVector(x)
robjects.globalEnv["nino3"] = nino3
#r script
r=robjects.r
rscripts="""
png('NINO3_R.png')
nino3ts=ts(nino3,start=c(%(year)s,%(month)s),frequency=12)
plot.ts(nino3ts)
dev.off()
"""
r(rscripts %locals())
rscript="""
new<-data.frame(rlevel=seq(-7,35,1)) # (1)
png(\'kuroshio_reg_R3.png\')
dc<-predict(fm,new,interval="confidence",level=0.99) #(2)
dp<-predict(fm,new,interval="prediction",level=0.99) #(3)
matplot(new$rlevel,cbind(dp,dc),lty=1,col=c("black","black","black","black","blue","blue"),type="l",xlab="",ylab="",xlim=c(-7, 35),ylim=c(30,33.5)) #(4)
par(new=T) #(5)
plot(rlevel,rlat,xlab=\"Kushimoto-Uragami\",ylab=\"Kuroshio Latitude\",xlim=c(-7, 35),ylim=c(30,33.5)) #(6)
dev.off()
"""
r(rscript)
図
r('png(\'kuroshio_reg_R2.png\')')
r('opar<-par(mfrow=c(2,2),oma=c(0,0,1.1,0))') #(1)
r('plot(fm)') #(2)
r('par(opar)') #(3)
r('dev.off()')
#read module
#for data
from read_kuroshio_lat import read_kuroshio_lat
from read_kushimoto_uragami import read_kushimoto_uragami
#for R
import rpy2.robjects as robjects
#read data
kuroshio_lat=read_kuroshio_lat()
kushimoto_uragami=read_kushimoto_uragami()
#for R
rlat=robjects.FloatVector(kuroshio_lat)
rlevel=robjects.FloatVector(kushimoto_uragami)
robjects.globalEnv["rlat"] = rlat
robjects.globalEnv["rlevel"] = rlevel
r=robjects.r
r('png(\'kuroshio_reg_R.png\')')
r('plot(rlevel,rlat,xlab=\"Kushimoto-Uragami\",ylab=\"Kuroshio Latitude\")') #scatter plot
r('fm<-lm(rlat ~ rlevel)') # linear model
r('abline(fm)') # plot regression line
r('dev.off()')
summary=r('sm<-summary(fm)')
print summary
図
print summary.names[1] "call" "terms" "residuals" "coefficients"
cf=r('cf<-sm$coefficients')[[1]]
print cf.names
for i in range(8):31.9559366577
print cf[i]
#boxplot
scripts="""
png(’boxplot.png’)
boxplot(rdatain)
dev.off()
"""
r(scripts)
#hist
scripts="""
png('hist.png')
hist(rdatain,seq(-4, 4, 0.4),prob=TRUE)
x<-seq(-4, 4, 0.2)
lines(x, dnorm(x, mean=mean(rdatain), sd=sqrt(var(rdatain))), lty=3)
dev.off()
"""
a=r(scripts)
#cdf plot
scripts="""
png('cdf.png')
plot(ecdf(rdatain), do.points=FALSE, verticals=TRUE)
x<-seq(-4, 4, 0.01)
lines(x, pnorm(x, mean=mean(rdatain), sd=sqrt(var(rdatain))), lty=3)
dev.off()
"""
r(scripts)
#qqplot
scripts="""
png('qqplot.png')
qqnorm(rdatain)
qqline(rdatain)
dev.off()
"""
r(scripts)
from readAO import readAO # read AO data
from pick_month import pick_month # pick month
#
import scikits.timeseries as ts
#
import rpy2.robjects as robjects
# read AO data
AOindex_series=readAO()
# pickup December
AOindex_dec=pick_month(AOindex_series,month=12)
#remove missig value
data=AOindex_dec[AOindex_dec.mask==False].data
# for R
rdata=robjects.FloatVector(data)
robjects.globalEnv["rdatain"] = rdata
r=robjects.r