Nothing
R/plot.etasclass.R
In etasFLP: Mixed FLP and ML Estimation of ETAS Space-Time Point Processes for Earthquake Description
Defines functions
plot.etasclass
Documented in
plot.etasclass
plot.etasclass
plot.etasclass
plot.etasclass <-
function(x,pdf=FALSE,file ="etasplot", ngrid=201,nclass=10,tfixed=0,flag.log=FALSE,...){
# function(x,pdf=FALSE,file ="etasplot", ngrid=201,flag.3D=FALSE,flag.log=FALSE,ellipse=FALSE,...){
if (!inherits(x,"etasclass"))stop("argument must be an etasclass object")
nsimps=1+(ngrid-1)/nclass
simpson=(nsimps==trunc(nsimps))
if(!simpson)stop("argument nclass must divide exactly ngrid-1 in order to use Simpson's rule to compute theoretical frequencies")
scaleres=TRUE
ellipse=FALSE
# kern.var=x$kern.var
kern.var=FALSE
# which=1:4
# maxgraph =7
# ind.graph =array(FALSE,maxgraph)
# ind.graph[which]=TRUE
#tfixed=0 add in the rd file
cat("Computation of the ETAS space intensity integrated intensity on a grid","\n")
cat("for large catalogs computation can take some minutes. Please wait ...","\n")
#####################################################
# plot time intensity # with points
typegraph=1
file=paste(file,".pdf",sep="")
if (pdf){
pdf(file=file,onefile=TRUE)}
#else
#{
#dev.new()
#}
#####################################################
### computation of back intensity and triggered intensity on a given x-y grid
### back.grid computed with the same kernel used for back.dens on observed points
## back.grid computed on km coordinates
# chamge the is.backconstant switch
magnitude=x$cat$magnitude
#if(!x$is.backconstant){
{
xcat.km =x$cat$xcat.work
ycat.km =x$cat$ycat.work
w =x$rho.weights
hdef =x$hdef
n =length(xcat.km)
tmax =x$tmax
params =x$params
mu= params[1]
k0= params[2]
c= params[3]
p= params[4]
# a= params[5]
gamma= params[5]
d= params[6]
q= params[7]
#
eps=1/n
rangex =range(xcat.km)
rangey =range(ycat.km)
eps.x =diff(rangex)*eps/2
eps.y =diff(rangey)*eps/2
rangex.km =c(min(rangex)-eps.x,max(rangex)+eps.x)
rangey.km =c(min(rangey)-eps.y,max(rangey)+eps.y)
ranget=diff(range(x$cat$time.work))
alpha =0
if (kern.var) alpha=hdef[3]
space.grid =xy.grid(rangex.km,rangey.km,ngrid)
ngridtot =length(space.grid[,1])
#kde=kde2dnew.fortran(xcat.km,ycat.km,space.grid[,1],space.grid[,2],w=w,factor.xy=1,h=hdef,kern.var=kern.var,alpha=alpha)
if(!x$is.backconstant){
kde=kde2dnew.fortran(xcat.km,ycat.km,space.grid[,1],space.grid[,2],w=w,factor.xy=1,h=hdef,hvarx=x$hvarx,hvary=x$hvary)
#wmat1 =matrix(kde$wmat,n,4)
back.grid =ranget*mu*kde$z
}
else
{
kde=kde2dnew.fortran(xcat.km,ycat.km,space.grid[,1],space.grid[,2],w=w,factor.xy=1,h=c(1,1))
#wmat1 =matrix(kde$wmat,n,4)
back.grid =ranget*mu*kde$z^0
}
#### computation MUST be in km IN THIS VERSION because parameters are evaluated in km
ris=.Fortran("etasfull8tintegratednew" ,NAOK=TRUE,
n=as.integer(n),
mu=as.double(mu),k=as.double(k0),
c=as.double(c),p=as.double(p),
g=as.double(gamma),
d=as.double(d),q=as.double(q),
x=as.double(xcat.km),y=as.double(ycat.km), t=as.double(x$cat$time.work),
m=as.double(magnitude),
predictor=as.double(x$predictor),
l=back.grid,
ngridtot=as.integer(ngridtot), xgrid=as.double(space.grid[,1]),
ygrid=as.double(space.grid[,2]),
tmax=as.double(tmax))
trig.grid =ris$l
### trig.grid intensity on a x-y grid by time integration of ETAS intensity function
###########################################################################################
}
tot.grid=back.grid+trig.grid
# maps of triggered intensity for a single day
totfixed.grid=tot.grid*0.
if(tfixed>0){
ind=x$cat$time.work<tfixed
ris=.Fortran("etasfull8tfixednew" ,NAOK=TRUE,
n=as.integer(sum(ind)),
mu=as.double(mu),k=as.double(k0),
c=as.double(c),p=as.double(p),
g=as.double(gamma),
d=as.double(d),q=as.double(q),
x=as.double(xcat.km[ind]),y=as.double(ycat.km[ind]), t=as.double(x$cat$time.work[ind]),m=as.double(magnitude[ind]),
predictor=as.double(x$predictor),
l=back.grid,
ngridtot=as.integer(ngridtot), xgrid=as.double(space.grid[,1]),
ygrid=as.double(space.grid[,2]),tfixed=as.double(tfixed))
totfixed.grid =ris$l+back.grid/ranget
if(flag.log)totfixed.grid=log(totfixed.grid)
}
if(flag.log) {
back.grid=log(back.grid)
trig.grid=log(trig.grid)
tot.grid=log( tot.grid)
}
## change grid to degrees
rangex =range(x$cat$long)
rangey =range(x$cat$lat)
eps.x =diff(rangex)*eps/2
eps.y =diff(rangey)*eps/2
rangex =c(min(rangex)-eps.x,max(rangex)+eps.x)
rangey =c(min(rangey)-eps.y,max(rangey)+eps.y)
space.grid =xy.grid(rangex,rangey,ngrid)
ngridtot =length(space.grid[,1])
x.grid =seq(rangex[1],rangex[2],length=ngrid)
y.grid =seq(rangey[1],rangey[2],length=ngrid)
if(tfixed>0){
### start triggered intensity plotting for a tfixed ########################################
if(!pdf) dev.new()
ind1=x$cat$time.work>=tfixed & x$cat$time.work<(tfixed+1)
mapxy=map("worldHires",xlim=range(x$cat$long),ylim=range(x$cat$lat),plot=FALSE)
image.plot(x.grid,y.grid,(matrix(totfixed.grid,ngrid,ngrid))
,col=gray.colors(128, start = 0., end = 1., gamma =2 )
,xlab="x-longitude",ylab="y-latitude"
,main=paste("Triggered intensity and observed points at day ",tfixed)
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green"
)
contour(x.grid,y.grid,(matrix(totfixed.grid,ngrid,ngrid)),col="red",add=TRUE)
points(x$cat$long[ind1],x$cat$lat[ind1],cex=sqrt(exp(x$cat$magnitude[ind1]))/8,col=4,pch=19)
### end triggered intensity plotting for a tfixed ########################################
}
if(!pdf) dev.new()
### start triggered intensity plotting
mapxy=map("worldHires",xlim=range(x$cat$long),ylim=range(x$cat$lat),plot=FALSE)
image.plot(x.grid,y.grid,(matrix(trig.grid,ngrid,ngrid))
,col=gray.colors(128, start = 0., end = 1., gamma =2 )
,xlab="x-longitude",ylab="y-latitude"
,main="Triggered Intensity"
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green"
)
contour(x.grid,y.grid,(matrix(trig.grid,ngrid,ngrid)),col="red",add=TRUE)
### end triggered intensity plotting
### start background intensity plotting
box()
if(!pdf) dev.new()
image.plot(x.grid,y.grid,(matrix(back.grid,ngrid,ngrid))
,col=gray.colors(128, start = 0., end = 1., gamma =2 ),
xlab="x-longitude",ylab="y-latitude",
main="Background Intensity"
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green")
contour(x.grid,y.grid,(matrix(back.grid,ngrid,ngrid)),col="red",add=TRUE)
box()
### start total intensity plotting
if(!pdf) dev.new()
image.plot(x.grid,y.grid,(matrix(tot.grid,ngrid,ngrid)),
col=gray.colors(128, start = 0., end = 1., gamma =2 ),
xlab="x-longitude",ylab="y-latitude",main="Total Intensity"
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green")
contour(x.grid,y.grid,(matrix(tot.grid,ngrid,ngrid)),col="red",add=TRUE)
box()
### start total intensity plotting with observed points
ts=(x$cat$time-min(x$cat$time))/diff(range(x$cat$time))
if(!pdf) dev.new()
image.plot(x.grid,y.grid,
(matrix(tot.grid,ngrid,ngrid)),
col=gray.colors(128, start = 0., end = 1., gamma =2 ),
xlab="x-longitude",ylab="y-latitude",main="Total Intensity with observed points \n Circles area proportional to magnitude; red: recent, blu:older"
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green")
points(x$cat$long,x$cat$lat,cex=sqrt(exp(x$cat$magnitude))/8,col=rgb(ts,0,1-ts),pch=19)
contour(x.grid,y.grid,(matrix(tot.grid,ngrid,ngrid)),col="yellow",add=TRUE)
box()
#####################################################
if((as.numeric(kern.var)*as.numeric(ellipse))==1){
wmat1 =matrix(0,n,4)
if(!pdf) dev.new()
wx=wmat1[,2]^2
wy=wmat1[,3]^2
wxy=wmat1[,4]*wmat1[,2]*wmat1[,3]
x=xcat.km
y=ycat.km
plot(x,y,type="p",pch=22,col=1,main="adaptive kernel")
for (j in 1:n)
cat(j)
# internal ellipses
lines(ellipse(matrix(c(wx[j],wxy[j],wxy[j],wy[j]),2,2),centre=c(x[j],y[j]),level=c(0.500)),col=2)
# external ellipses
for (j in 1:n) lines(ellipse(matrix(c(wx[j],wxy[j],wxy[j],wy[j]),2,2),centre=c(x[j],y[j]),level=c(0.999)),col=3)
}
etas.l=x$l
if(flag.log) etas.l=log(etas.l)
#if(flag.3D){
#typegraph=2
#plot3d(x$cat$long,x$cat$lat,etas.l,type="n",zlab=paste("estimated intensity ","lambda(x,y)"),
#xlab="x-longitude",ylab="y-latitude",main="Estimated intensities in observed points")
#lines3d(cbind(mapxy$x,mapxy$y,min(etas.l)),col="red")
#plot3d(x$cat$long,x$cat$lat,etas.l,add=TRUE, type="h")
#}
#####################################################################################
### added november 7th, 2014
# computation of intensities in observed points integrated with respect to time
#
ngridtot =n
back.obs =ranget*mu*kde2dnew.fortran(xcat.km,ycat.km,xcat.km,ycat.km,w=w,factor.xy=1,h=hdef,hvarx=x$hvarx,hvary=x$hvary)$z
#### computation MUST be in km IN THIS VERSION because parameters are evaluated in km
ris=.Fortran("etasfull8tintegratednew" ,NAOK=TRUE,
n=as.integer(n),
mu=as.double(mu),k=as.double(k0),
c=as.double(c),p=as.double(p),
g=as.double(gamma),
d=as.double(d),q=as.double(q),
x=as.double(xcat.km),y=as.double(ycat.km), t=as.double(x$cat$time.work),m=as.double(x$cat$magn1),
predictor=as.double(x$predictor),
l=back.obs,
ngridtot=as.integer(n), xgrid=as.double(xcat.km), ygrid=as.double(ycat.km),tmax=as.double(tmax))
trig.obs =ris$l
tot.obs=back.obs+trig.obs
teo1=matrix(0,nclass,nclass)
teo2=matrix(0,nclass,nclass)
teo3=matrix(0,nclass,nclass)
teo4=matrix(0,nclass,nclass)
emp1=teo1
emp2=teo2
emp3=teo3
emp4=teo4
coeff.integr =simpson.kD(nsimps,2)
dx = diff(rangex.km)/(ngrid-1)
dy = diff(rangey.km)/(ngrid-1)
coeff.integr = as.vector(coeff.integr*dx*dy)
x0=min(rangex.km)
y0=min(rangey.km)
xint=diff(rangex.km)/nclass
yint=diff(rangey.km)/nclass
for (i in 1:nclass){
for (j in 1:nclass){
ivert1=(nsimps-1)*(i-1)+1
ivert2=(nsimps-1)*i+1
jvert1=(nsimps-1)*(j-1)+1
jvert2=(nsimps-1)*j+1
xymat=expand.grid(ivert1:ivert2,jvert1:jvert2)
xyvec=ngrid*(xymat[,2]-1)+xymat[,1]
teovec1=tot.grid[xyvec]
teovec2=back.grid[xyvec]
teo1[i,j]=sum(coeff.integr*teovec1)
#teo2[i,j]=sum(coeff.integr*sqrt(teovec))
teo2[i,j]=sum(coeff.integr*teovec2)
x1=x0+(i-1)*xint
x2=x1+xint
y1=y0+(j-1)*yint
y2=y1+yint
indrid=(xcat.km >= x1)&(xcat.km < x2)&(ycat.km >= y1)&(ycat.km < y2)
teo3[i,j]=sum(coeff.integr*sqrt(teovec1))
teo4[i,j]=sum(coeff.integr)
emp1[i,j]=sum(indrid)
emp2[i,j]=sum(x$rho.weights[indrid])
emp3[i,j]=sum(1/sqrt(tot.obs[indrid]))
emp4[i,j]=0
if(emp1[i,j]>0) emp4[i,j]=sum(1/tot.obs[indrid])
}
}
#####################################################################################
## plotting space residuals
xx=range(x.grid)
yy=range(y.grid)
hx=diff(xx)/nclass
hy=diff(yy)/nclass
x.grid =seq(xx[1]+hx/2,xx[2]-hx/2,length=nclass)
y.grid =seq(yy[1]+hy/2,yy[2]-hy/2,length=nclass)
alpha=255
cblue= rgb(255:0, 255:0, 255, maxColorValue = 255,alpha=alpha)
cred = rgb(255, 255:0, 255:0, maxColorValue = 255,alpha=alpha)
bluered=c(cred[256:1],cblue)
smax=3
std=(emp1-teo1)/sqrt(teo1)
std4=std<(-smax)
std[std4]=-smax
std4=std>(smax)
std[std4]=smax
if(!pdf) dev.new()
image.plot(x.grid,y.grid,std,zlim=c(-smax,smax),col=bluered,
main="Standardized differences between \n observed and theoretical frequency (whole model)")
#image.plot(x.grid,y.grid,std,zlim=c(-smax,smax),col=bluered,xaxp=c(xx[1],xx[2],nclass),yaxp=c(yy[1],yy[2],nclass))
grid(nclass)
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",col="black")
############### for background
std=(emp2-teo2)/sqrt(teo2)
std4=std<(-smax)
std[std4]=-smax
std4=std>(smax)
std[std4]=smax
if(!pdf) dev.new()
image.plot(x.grid,y.grid,std,zlim=c(-smax,smax),col=bluered,
main="Standardized differences between \n observed and theoretical frequency (background only)")
#image.plot(x.grid,y.grid,std,zlim=c(-smax,smax),col=bluered,xaxp=c(xx[1],xx[2],nclass),yaxp=c(yy[1],yy[2],nclass))
grid(nclass)
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",col="black")
############################
######### four graph on a page
if(!pdf) dev.new()
par(mfcol=c(3,2))
## plotting observed vs. observed theoretical
np=200
# page7: graph2
par(cex.main=0.8,mar=c(6, 4, 5, 2) + 0.2)
plot(teo1,emp1,xlab="theor. freq",ylab="empir. freq",main="Total space intens.")
abline(a=0,b=1)
tvec=seq(0,max(teo1),length.out=np)
yvec=1.96*sqrt(tvec)
lines(tvec,tvec+yvec,col=2)
lines(tvec,tvec-yvec,col=2)
##
## plotting space residuals
xx=range(x.grid)
yy=range(y.grid)
hx=diff(xx)/nclass
hy=diff(yy)/nclass
# page7: graph3
plot(teo2,emp2,xlab="theor. freq",ylab="empir. freq",main="Backgr. space intens.")
abline(a=0,b=1)
tvec=seq(0,max(teo2),length.out=np)
yvec=1.96*sqrt(tvec)
lines(tvec,tvec+yvec,col=2)
lines(tvec,tvec-yvec,col=2)
# page7: graph4
plot(teo1-teo2,emp1-emp2,,xlab="theor. freq",ylab="empir. freq",main="Triggered space intens.")
abline(a=0,b=1)
tvec=seq(0,max(teo1-teo2),length.out=np)
yvec=1.96*sqrt(tvec)
lines(tvec,tvec+yvec,col=2)
lines(tvec,tvec-yvec,col=2)
### page 7 scaled residuals (second column of graph)
# page8: graph1
## qqplot
#std=(emp1-teo1)/sqrt(teo1)
#qqnorm(std); qqline(std, col = 2)
## plotting residuals vs. theoretical
## on transformed scales
std=(emp1-teo1)/sqrt(teo1)
# page7: graph1
plot(sqrt(teo1),std,xlab="sqrt theor. freq",ylab="std.zd residuals",main="Total space intens.")
abline(a=-1.96,b=0,col=2)
abline(a=+1.96,b=0,col=2)
##
## plotting space residuals
## plotting residuals vs. theoretical
## on transformed scales
std=(emp2-teo2)/sqrt(teo2)
# page7: graph3
plot(sqrt(teo2),std,xlab="sqrt theor. freq",ylab="std.zd residuals",main="Background space intens.")
abline(a=-1.96,b=0,col=2)
abline(a=+1.96,b=0,col=2)
## plotting residuals vs. theoretical
## on transformed scales
std=(emp1-emp2-(teo1-teo2))/sqrt(teo1-teo2)
# page8: graph4
plot(sqrt(teo1-teo2),std,xlab="sqrt theor. freq",ylab="std.zd residuals",main="Triggered space intens.")
abline(a=-1.96,b=0,col=2)
abline(a=+1.96,b=0,col=2)
############################## time residuals #######################################
######## 9-12-2014 #############
## transformed points. intensity integrated with respect to space ####
cat("Computation of the time intensity transform","\n")
lambda=mu
ntemp =n
iprint =FALSE
ttot =x$cat$time.work
rho.tot =x$rho.s2
t.trasf =array(0,ntemp)
t.intens=array(0,ntemp)
t.intens[1]=x$back.dens[1]*lambda
ntheta =ncol(rho.tot)
for (i in 2:ntemp){
times =ttot[1:i]
range.t =diff(range(times))
rho.s2 =rho.tot[1:i,]
tmax.i =max(times)
magnitudes=magnitude[1:i]
predictor.i=x$predictor[1:i]
trig=(c+tmax.i-times)^(-p)
if(p==1){
it=log(c+tmax.i-times)-log(c)
}
else
{
it=((c+tmax.i-times)^(1-p)-c^(1-p))/(1-p)
}
if(x$onlytime){
integral=k0*sum(exp(predictor.i)*it)
}
else
##############################################################################
# starting space integration (polar transformation)
##############################################################################
{
m1 =as.vector(exp(predictor.i))
m2 =as.vector(exp(gamma*magnitudes))
### approximate polar integration to whole space
### approximate polar integration to whole space by division in ntheta triangles centered in xi,yi
etas =rowSums(m1*m2*((rho.s2*rho.s2/m2+d)^(1-q)-d^(1-q)))
space =(pi/((1-q)*ntheta))*etas
integral=sum(k0*it*space)
integral.intens=sum(k0*trig*space)
integraltot =integral +lambda*range.t*x$back.integral
integraltot.intens =integral.intens+lambda*x$back.dens[i]
if(iprint) cat(i," ")
t.trasf[i] =integraltot
t.intens[i] =integraltot.intens
}
##############################################################################
#
# end space integration
#
##############################################################################
}
if(!pdf) dev.new()
#close.screen()
par(mfcol=c(2,1))
#split.screen(rbind(c(0,1,0.4,1),c(0,1,0,0.4)))
#screen(1)
plot(t.trasf,type="l",xlab="i",ylab=expression(tau(t[i])),
main="Transformed times")
abline(a=0,b=1,col=2)
#screen(2)
plot(t.intens,type="h",xlab=expression(t[i]),ylab=expression(lambda(t[i])))
close.screen()
################################ end of time residuals
########## end of graphic output ################
if(pdf) dev.off()
return(list(
x.grid=x.grid,
y.grid=y.grid,
space.grid=space.grid,
back.grid=back.grid,
trig.grid=trig.grid,
tot.grid=tot.grid,
tfixed=tfixed,
totfixed.grid=totfixed.grid,
back.obs=back.obs,
trig.obs=trig.obs,
tot.obs=tot.obs,
teo1=teo1,
teo2=teo2,
#teo3=teo3,
#teo4=teo4,
emp1=emp1,
emp2=emp2,
#emp3=emp3,
#emp4=emp4,
dx=dx,
dy=dy,
ranget=ranget,
# etasn=etasn,
magnitude=magnitude,
t.trasf=t.trasf
#,wmat
))
}
Try the etasFLP package in your browser
Any scripts or data that you put into this service are public.
etasFLP documentation built on Sept. 14, 2023, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot.etasclass
Documented in plot.etasclass plot.etasclass plot.etasclass
plot.etasclass <-
function(x,pdf=FALSE,file ="etasplot", ngrid=201,nclass=10,tfixed=0,flag.log=FALSE,...){
# function(x,pdf=FALSE,file ="etasplot", ngrid=201,flag.3D=FALSE,flag.log=FALSE,ellipse=FALSE,...){
if (!inherits(x,"etasclass"))stop("argument must be an etasclass object")
nsimps=1+(ngrid-1)/nclass
simpson=(nsimps==trunc(nsimps))
if(!simpson)stop("argument nclass must divide exactly ngrid-1 in order to use Simpson's rule to compute theoretical frequencies")
scaleres=TRUE
ellipse=FALSE
# kern.var=x$kern.var
kern.var=FALSE
# which=1:4
# maxgraph =7
# ind.graph =array(FALSE,maxgraph)
# ind.graph[which]=TRUE
#tfixed=0 add in the rd file
cat("Computation of the ETAS space intensity integrated intensity on a grid","\n")
cat("for large catalogs computation can take some minutes. Please wait ...","\n")
#####################################################
# plot time intensity # with points
typegraph=1
file=paste(file,".pdf",sep="")
if (pdf){
pdf(file=file,onefile=TRUE)}
#else
#{
#dev.new()
#}
#####################################################
### computation of back intensity and triggered intensity on a given x-y grid
### back.grid computed with the same kernel used for back.dens on observed points
## back.grid computed on km coordinates
# chamge the is.backconstant switch
magnitude=x$cat$magnitude
#if(!x$is.backconstant){
{
xcat.km =x$cat$xcat.work
ycat.km =x$cat$ycat.work
w =x$rho.weights
hdef =x$hdef
n =length(xcat.km)
tmax =x$tmax
params =x$params
mu= params[1]
k0= params[2]
c= params[3]
p= params[4]
# a= params[5]
gamma= params[5]
d= params[6]
q= params[7]
#
eps=1/n
rangex =range(xcat.km)
rangey =range(ycat.km)
eps.x =diff(rangex)*eps/2
eps.y =diff(rangey)*eps/2
rangex.km =c(min(rangex)-eps.x,max(rangex)+eps.x)
rangey.km =c(min(rangey)-eps.y,max(rangey)+eps.y)
ranget=diff(range(x$cat$time.work))
alpha =0
if (kern.var) alpha=hdef[3]
space.grid =xy.grid(rangex.km,rangey.km,ngrid)
ngridtot =length(space.grid[,1])
#kde=kde2dnew.fortran(xcat.km,ycat.km,space.grid[,1],space.grid[,2],w=w,factor.xy=1,h=hdef,kern.var=kern.var,alpha=alpha)
if(!x$is.backconstant){
kde=kde2dnew.fortran(xcat.km,ycat.km,space.grid[,1],space.grid[,2],w=w,factor.xy=1,h=hdef,hvarx=x$hvarx,hvary=x$hvary)
#wmat1 =matrix(kde$wmat,n,4)
back.grid =ranget*mu*kde$z
}
else
{
kde=kde2dnew.fortran(xcat.km,ycat.km,space.grid[,1],space.grid[,2],w=w,factor.xy=1,h=c(1,1))
#wmat1 =matrix(kde$wmat,n,4)
back.grid =ranget*mu*kde$z^0
}
#### computation MUST be in km IN THIS VERSION because parameters are evaluated in km
ris=.Fortran("etasfull8tintegratednew" ,NAOK=TRUE,
n=as.integer(n),
mu=as.double(mu),k=as.double(k0),
c=as.double(c),p=as.double(p),
g=as.double(gamma),
d=as.double(d),q=as.double(q),
x=as.double(xcat.km),y=as.double(ycat.km), t=as.double(x$cat$time.work),
m=as.double(magnitude),
predictor=as.double(x$predictor),
l=back.grid,
ngridtot=as.integer(ngridtot), xgrid=as.double(space.grid[,1]),
ygrid=as.double(space.grid[,2]),
tmax=as.double(tmax))
trig.grid =ris$l
### trig.grid intensity on a x-y grid by time integration of ETAS intensity function
###########################################################################################
}
tot.grid=back.grid+trig.grid
# maps of triggered intensity for a single day
totfixed.grid=tot.grid*0.
if(tfixed>0){
ind=x$cat$time.work<tfixed
ris=.Fortran("etasfull8tfixednew" ,NAOK=TRUE,
n=as.integer(sum(ind)),
mu=as.double(mu),k=as.double(k0),
c=as.double(c),p=as.double(p),
g=as.double(gamma),
d=as.double(d),q=as.double(q),
x=as.double(xcat.km[ind]),y=as.double(ycat.km[ind]), t=as.double(x$cat$time.work[ind]),m=as.double(magnitude[ind]),
predictor=as.double(x$predictor),
l=back.grid,
ngridtot=as.integer(ngridtot), xgrid=as.double(space.grid[,1]),
ygrid=as.double(space.grid[,2]),tfixed=as.double(tfixed))
totfixed.grid =ris$l+back.grid/ranget
if(flag.log)totfixed.grid=log(totfixed.grid)
}
if(flag.log) {
back.grid=log(back.grid)
trig.grid=log(trig.grid)
tot.grid=log( tot.grid)
}
## change grid to degrees
rangex =range(x$cat$long)
rangey =range(x$cat$lat)
eps.x =diff(rangex)*eps/2
eps.y =diff(rangey)*eps/2
rangex =c(min(rangex)-eps.x,max(rangex)+eps.x)
rangey =c(min(rangey)-eps.y,max(rangey)+eps.y)
space.grid =xy.grid(rangex,rangey,ngrid)
ngridtot =length(space.grid[,1])
x.grid =seq(rangex[1],rangex[2],length=ngrid)
y.grid =seq(rangey[1],rangey[2],length=ngrid)
if(tfixed>0){
### start triggered intensity plotting for a tfixed ########################################
if(!pdf) dev.new()
ind1=x$cat$time.work>=tfixed & x$cat$time.work<(tfixed+1)
mapxy=map("worldHires",xlim=range(x$cat$long),ylim=range(x$cat$lat),plot=FALSE)
image.plot(x.grid,y.grid,(matrix(totfixed.grid,ngrid,ngrid))
,col=gray.colors(128, start = 0., end = 1., gamma =2 )
,xlab="x-longitude",ylab="y-latitude"
,main=paste("Triggered intensity and observed points at day ",tfixed)
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green"
)
contour(x.grid,y.grid,(matrix(totfixed.grid,ngrid,ngrid)),col="red",add=TRUE)
points(x$cat$long[ind1],x$cat$lat[ind1],cex=sqrt(exp(x$cat$magnitude[ind1]))/8,col=4,pch=19)
### end triggered intensity plotting for a tfixed ########################################
}
if(!pdf) dev.new()
### start triggered intensity plotting
mapxy=map("worldHires",xlim=range(x$cat$long),ylim=range(x$cat$lat),plot=FALSE)
image.plot(x.grid,y.grid,(matrix(trig.grid,ngrid,ngrid))
,col=gray.colors(128, start = 0., end = 1., gamma =2 )
,xlab="x-longitude",ylab="y-latitude"
,main="Triggered Intensity"
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green"
)
contour(x.grid,y.grid,(matrix(trig.grid,ngrid,ngrid)),col="red",add=TRUE)
### end triggered intensity plotting
### start background intensity plotting
box()
if(!pdf) dev.new()
image.plot(x.grid,y.grid,(matrix(back.grid,ngrid,ngrid))
,col=gray.colors(128, start = 0., end = 1., gamma =2 ),
xlab="x-longitude",ylab="y-latitude",
main="Background Intensity"
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green")
contour(x.grid,y.grid,(matrix(back.grid,ngrid,ngrid)),col="red",add=TRUE)
box()
### start total intensity plotting
if(!pdf) dev.new()
image.plot(x.grid,y.grid,(matrix(tot.grid,ngrid,ngrid)),
col=gray.colors(128, start = 0., end = 1., gamma =2 ),
xlab="x-longitude",ylab="y-latitude",main="Total Intensity"
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green")
contour(x.grid,y.grid,(matrix(tot.grid,ngrid,ngrid)),col="red",add=TRUE)
box()
### start total intensity plotting with observed points
ts=(x$cat$time-min(x$cat$time))/diff(range(x$cat$time))
if(!pdf) dev.new()
image.plot(x.grid,y.grid,
(matrix(tot.grid,ngrid,ngrid)),
col=gray.colors(128, start = 0., end = 1., gamma =2 ),
xlab="x-longitude",ylab="y-latitude",main="Total Intensity with observed points \n Circles area proportional to magnitude; red: recent, blu:older"
)
grid(col="grey")
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",
xlim=range(x$cat$long),ylim=range(x$cat$lat),col="green")
points(x$cat$long,x$cat$lat,cex=sqrt(exp(x$cat$magnitude))/8,col=rgb(ts,0,1-ts),pch=19)
contour(x.grid,y.grid,(matrix(tot.grid,ngrid,ngrid)),col="yellow",add=TRUE)
box()
#####################################################
if((as.numeric(kern.var)*as.numeric(ellipse))==1){
wmat1 =matrix(0,n,4)
if(!pdf) dev.new()
wx=wmat1[,2]^2
wy=wmat1[,3]^2
wxy=wmat1[,4]*wmat1[,2]*wmat1[,3]
x=xcat.km
y=ycat.km
plot(x,y,type="p",pch=22,col=1,main="adaptive kernel")
for (j in 1:n)
cat(j)
# internal ellipses
lines(ellipse(matrix(c(wx[j],wxy[j],wxy[j],wy[j]),2,2),centre=c(x[j],y[j]),level=c(0.500)),col=2)
# external ellipses
for (j in 1:n) lines(ellipse(matrix(c(wx[j],wxy[j],wxy[j],wy[j]),2,2),centre=c(x[j],y[j]),level=c(0.999)),col=3)
}
etas.l=x$l
if(flag.log) etas.l=log(etas.l)
#if(flag.3D){
#typegraph=2
#plot3d(x$cat$long,x$cat$lat,etas.l,type="n",zlab=paste("estimated intensity ","lambda(x,y)"),
#xlab="x-longitude",ylab="y-latitude",main="Estimated intensities in observed points")
#lines3d(cbind(mapxy$x,mapxy$y,min(etas.l)),col="red")
#plot3d(x$cat$long,x$cat$lat,etas.l,add=TRUE, type="h")
#}
#####################################################################################
### added november 7th, 2014
# computation of intensities in observed points integrated with respect to time
#
ngridtot =n
back.obs =ranget*mu*kde2dnew.fortran(xcat.km,ycat.km,xcat.km,ycat.km,w=w,factor.xy=1,h=hdef,hvarx=x$hvarx,hvary=x$hvary)$z
#### computation MUST be in km IN THIS VERSION because parameters are evaluated in km
ris=.Fortran("etasfull8tintegratednew" ,NAOK=TRUE,
n=as.integer(n),
mu=as.double(mu),k=as.double(k0),
c=as.double(c),p=as.double(p),
g=as.double(gamma),
d=as.double(d),q=as.double(q),
x=as.double(xcat.km),y=as.double(ycat.km), t=as.double(x$cat$time.work),m=as.double(x$cat$magn1),
predictor=as.double(x$predictor),
l=back.obs,
ngridtot=as.integer(n), xgrid=as.double(xcat.km), ygrid=as.double(ycat.km),tmax=as.double(tmax))
trig.obs =ris$l
tot.obs=back.obs+trig.obs
teo1=matrix(0,nclass,nclass)
teo2=matrix(0,nclass,nclass)
teo3=matrix(0,nclass,nclass)
teo4=matrix(0,nclass,nclass)
emp1=teo1
emp2=teo2
emp3=teo3
emp4=teo4
coeff.integr =simpson.kD(nsimps,2)
dx = diff(rangex.km)/(ngrid-1)
dy = diff(rangey.km)/(ngrid-1)
coeff.integr = as.vector(coeff.integr*dx*dy)
x0=min(rangex.km)
y0=min(rangey.km)
xint=diff(rangex.km)/nclass
yint=diff(rangey.km)/nclass
for (i in 1:nclass){
for (j in 1:nclass){
ivert1=(nsimps-1)*(i-1)+1
ivert2=(nsimps-1)*i+1
jvert1=(nsimps-1)*(j-1)+1
jvert2=(nsimps-1)*j+1
xymat=expand.grid(ivert1:ivert2,jvert1:jvert2)
xyvec=ngrid*(xymat[,2]-1)+xymat[,1]
teovec1=tot.grid[xyvec]
teovec2=back.grid[xyvec]
teo1[i,j]=sum(coeff.integr*teovec1)
#teo2[i,j]=sum(coeff.integr*sqrt(teovec))
teo2[i,j]=sum(coeff.integr*teovec2)
x1=x0+(i-1)*xint
x2=x1+xint
y1=y0+(j-1)*yint
y2=y1+yint
indrid=(xcat.km >= x1)&(xcat.km < x2)&(ycat.km >= y1)&(ycat.km < y2)
teo3[i,j]=sum(coeff.integr*sqrt(teovec1))
teo4[i,j]=sum(coeff.integr)
emp1[i,j]=sum(indrid)
emp2[i,j]=sum(x$rho.weights[indrid])
emp3[i,j]=sum(1/sqrt(tot.obs[indrid]))
emp4[i,j]=0
if(emp1[i,j]>0) emp4[i,j]=sum(1/tot.obs[indrid])
}
}
#####################################################################################
## plotting space residuals
xx=range(x.grid)
yy=range(y.grid)
hx=diff(xx)/nclass
hy=diff(yy)/nclass
x.grid =seq(xx[1]+hx/2,xx[2]-hx/2,length=nclass)
y.grid =seq(yy[1]+hy/2,yy[2]-hy/2,length=nclass)
alpha=255
cblue= rgb(255:0, 255:0, 255, maxColorValue = 255,alpha=alpha)
cred = rgb(255, 255:0, 255:0, maxColorValue = 255,alpha=alpha)
bluered=c(cred[256:1],cblue)
smax=3
std=(emp1-teo1)/sqrt(teo1)
std4=std<(-smax)
std[std4]=-smax
std4=std>(smax)
std[std4]=smax
if(!pdf) dev.new()
image.plot(x.grid,y.grid,std,zlim=c(-smax,smax),col=bluered,
main="Standardized differences between \n observed and theoretical frequency (whole model)")
#image.plot(x.grid,y.grid,std,zlim=c(-smax,smax),col=bluered,xaxp=c(xx[1],xx[2],nclass),yaxp=c(yy[1],yy[2],nclass))
grid(nclass)
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",col="black")
############### for background
std=(emp2-teo2)/sqrt(teo2)
std4=std<(-smax)
std[std4]=-smax
std4=std>(smax)
std[std4]=smax
if(!pdf) dev.new()
image.plot(x.grid,y.grid,std,zlim=c(-smax,smax),col=bluered,
main="Standardized differences between \n observed and theoretical frequency (background only)")
#image.plot(x.grid,y.grid,std,zlim=c(-smax,smax),col=bluered,xaxp=c(xx[1],xx[2],nclass),yaxp=c(yy[1],yy[2],nclass))
grid(nclass)
map("worldHires",add=TRUE,xlab="x-longitude",ylab="y-latitude",col="black")
############################
######### four graph on a page
if(!pdf) dev.new()
par(mfcol=c(3,2))
## plotting observed vs. observed theoretical
np=200
# page7: graph2
par(cex.main=0.8,mar=c(6, 4, 5, 2) + 0.2)
plot(teo1,emp1,xlab="theor. freq",ylab="empir. freq",main="Total space intens.")
abline(a=0,b=1)
tvec=seq(0,max(teo1),length.out=np)
yvec=1.96*sqrt(tvec)
lines(tvec,tvec+yvec,col=2)
lines(tvec,tvec-yvec,col=2)
##
## plotting space residuals
xx=range(x.grid)
yy=range(y.grid)
hx=diff(xx)/nclass
hy=diff(yy)/nclass
# page7: graph3
plot(teo2,emp2,xlab="theor. freq",ylab="empir. freq",main="Backgr. space intens.")
abline(a=0,b=1)
tvec=seq(0,max(teo2),length.out=np)
yvec=1.96*sqrt(tvec)
lines(tvec,tvec+yvec,col=2)
lines(tvec,tvec-yvec,col=2)
# page7: graph4
plot(teo1-teo2,emp1-emp2,,xlab="theor. freq",ylab="empir. freq",main="Triggered space intens.")
abline(a=0,b=1)
tvec=seq(0,max(teo1-teo2),length.out=np)
yvec=1.96*sqrt(tvec)
lines(tvec,tvec+yvec,col=2)
lines(tvec,tvec-yvec,col=2)
### page 7 scaled residuals (second column of graph)
# page8: graph1
## qqplot
#std=(emp1-teo1)/sqrt(teo1)
#qqnorm(std); qqline(std, col = 2)
## plotting residuals vs. theoretical
## on transformed scales
std=(emp1-teo1)/sqrt(teo1)
# page7: graph1
plot(sqrt(teo1),std,xlab="sqrt theor. freq",ylab="std.zd residuals",main="Total space intens.")
abline(a=-1.96,b=0,col=2)
abline(a=+1.96,b=0,col=2)
##
## plotting space residuals
## plotting residuals vs. theoretical
## on transformed scales
std=(emp2-teo2)/sqrt(teo2)
# page7: graph3
plot(sqrt(teo2),std,xlab="sqrt theor. freq",ylab="std.zd residuals",main="Background space intens.")
abline(a=-1.96,b=0,col=2)
abline(a=+1.96,b=0,col=2)
## plotting residuals vs. theoretical
## on transformed scales
std=(emp1-emp2-(teo1-teo2))/sqrt(teo1-teo2)
# page8: graph4
plot(sqrt(teo1-teo2),std,xlab="sqrt theor. freq",ylab="std.zd residuals",main="Triggered space intens.")
abline(a=-1.96,b=0,col=2)
abline(a=+1.96,b=0,col=2)
############################## time residuals #######################################
######## 9-12-2014 #############
## transformed points. intensity integrated with respect to space ####
cat("Computation of the time intensity transform","\n")
lambda=mu
ntemp =n
iprint =FALSE
ttot =x$cat$time.work
rho.tot =x$rho.s2
t.trasf =array(0,ntemp)
t.intens=array(0,ntemp)
t.intens[1]=x$back.dens[1]*lambda
ntheta =ncol(rho.tot)
for (i in 2:ntemp){
times =ttot[1:i]
range.t =diff(range(times))
rho.s2 =rho.tot[1:i,]
tmax.i =max(times)
magnitudes=magnitude[1:i]
predictor.i=x$predictor[1:i]
trig=(c+tmax.i-times)^(-p)
if(p==1){
it=log(c+tmax.i-times)-log(c)
}
else
{
it=((c+tmax.i-times)^(1-p)-c^(1-p))/(1-p)
}
if(x$onlytime){
integral=k0*sum(exp(predictor.i)*it)
}
else
##############################################################################
# starting space integration (polar transformation)
##############################################################################
{
m1 =as.vector(exp(predictor.i))
m2 =as.vector(exp(gamma*magnitudes))
### approximate polar integration to whole space
### approximate polar integration to whole space by division in ntheta triangles centered in xi,yi
etas =rowSums(m1*m2*((rho.s2*rho.s2/m2+d)^(1-q)-d^(1-q)))
space =(pi/((1-q)*ntheta))*etas
integral=sum(k0*it*space)
integral.intens=sum(k0*trig*space)
integraltot =integral +lambda*range.t*x$back.integral
integraltot.intens =integral.intens+lambda*x$back.dens[i]
if(iprint) cat(i," ")
t.trasf[i] =integraltot
t.intens[i] =integraltot.intens
}
##############################################################################
#
# end space integration
#
##############################################################################
}
if(!pdf) dev.new()
#close.screen()
par(mfcol=c(2,1))
#split.screen(rbind(c(0,1,0.4,1),c(0,1,0,0.4)))
#screen(1)
plot(t.trasf,type="l",xlab="i",ylab=expression(tau(t[i])),
main="Transformed times")
abline(a=0,b=1,col=2)
#screen(2)
plot(t.intens,type="h",xlab=expression(t[i]),ylab=expression(lambda(t[i])))
close.screen()
################################ end of time residuals
########## end of graphic output ################
if(pdf) dev.off()
return(list(
x.grid=x.grid,
y.grid=y.grid,
space.grid=space.grid,
back.grid=back.grid,
trig.grid=trig.grid,
tot.grid=tot.grid,
tfixed=tfixed,
totfixed.grid=totfixed.grid,
back.obs=back.obs,
trig.obs=trig.obs,
tot.obs=tot.obs,
teo1=teo1,
teo2=teo2,
#teo3=teo3,
#teo4=teo4,
emp1=emp1,
emp2=emp2,
#emp3=emp3,
#emp4=emp4,
dx=dx,
dy=dy,
ranget=ranget,
# etasn=etasn,
magnitude=magnitude,
t.trasf=t.trasf
#,wmat
))
}
Try the etasFLP package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.