n=10000
QMAP=matrix(0,n,3)
MAPmat=matrix(0,n,N)
for(j in 1:N){
MAPmat[,j]=matrix(MAPS[[j]],n,1)}
for(i in 1:n){QMAP[i,]=createQuantiles(MAPmap[i,])}
QMAPf=NULL
for(i in 1:3){
QMAPf[[i]]=matrix(QMAP[,i],100,100)}
return(QMAPf)}
createQuantiles=function(vect,flag=0){
a=numeric(3)
a[1]=quantile(vect,0.05)
if(!flag){a[2]=mean(vect)} 
if(flag){a[2]=quantile(vect,0.5)}
a[3]=quantile(vect,0.95)
return(a)}
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
for(i in 1:5){a=createQuantiles(coeff[,i])}
return(a)}
 MonteCarloCoeff(EXP[[1]],10^4)
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
a=numeric(5)
for(i in 1:5){a[i]=createQuantiles(coeff[,i])}
return(a)}
 MonteCarloCoeff(EXP[[1]],10^4)
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
a=matrix(0,3,5)
for(i in 1:5){a[,i]=createQuantiles(coeff[,i])}
return(a)}
MonteCarloCoeff(EXP[[1]],10^4)
?matrix()
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
coln=c('Homog','EpoI_II','EpoIII','Faults,'Fracts'); rown=c('5%ile',mean,'95%ile')
a=matrix(0,3,5,dimnamnes=list(coln,rown))
for(i in 1:5){a[,i]=createQuantiles(coeff[,i])}
return(a)}
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
coln=c('Homog','EpoI_II','EpoIII','Faults,'Fracts'); rown=c('5%ile','mean','95%ile')
a=matrix(0,3,5,dimnamnes=list(coln,rown))
for(i in 1:5){a[,i]=createQuantiles(coeff[,i])}
return(a)}
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
coln=c('Homog','EpoI_II','EpoIII','Faults','Fracts'); rown=c('5%ile','mean','95%ile')
a=matrix(0,3,5,dimnamnes=list(coln,rown))
for(i in 1:5){a[,i]=createQuantiles(coeff[,i])}
return(a)}
return(a)}
MonteCarloCoeff(EXP[[1]],10^4)
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
coln=c('Homog','EpoI_II','EpoIII','Faults','Fracts'); rown=c('5%ile','mean','95%ile')
a=matrix(0,3,5,dimnames=list(coln,rown))
for(i in 1:5){a[,i]=createQuantiles(coeff[,i])}
return(a)}
MonteCarloCoeff(EXP[[1]],10^4)
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
coln=c('Homog','EpoI_II','EpoIII','Faults','Fracts'); rown=c('5%ile','mean','95%ile')
a=matrix(0,3,5,dimnames=list(rown,coln))
for(i in 1:5){a[,i]=createQuantiles(coeff[,i])}
return(a)}
MonteCarloCoeff(EXP[[1]],10^4)
seekMAPS=function(MAPS,N){ 
n=length(MAPS[[1]])
QMAP=matrix(0,n,3)
MAPmat=matrix(0,n,N)
for(j in 1:N){
MAPmat[,j]=matrix(MAPS[[j]],n,1)}
for(i in 1:n){QMAP[i,]=createQuantiles(MAPmap[i,])}
QMAPf=NULL
for(i in 1:3){
QMAPf[[i]]=matrix(QMAP[,i],100,100)}
return(QMAPf)}
MonteCarloMap(EXP[[1]],10^4)
seekMAPS=function(MAPS,N){ 
n=length(MAPS[[1]])
QMAP=matrix(0,n,3)
MAPmat=matrix(0,n,N)
for(j in 1:N){
MAPmat[,j]=matrix(MAPS[[j]],n,1)}
for(i in 1:n){QMAP[i,]=createQuantiles(MAPmat[i,])}
QMAPf=NULL
for(i in 1:3){
QMAPf[[i]]=matrix(QMAP[,i],100,100)}
return(QMAPf)}
MonteCarloMap(EXP[[1]],10^4)
outp=MonteCarloMap(EXP[[1]],10^4)
MonteCarloMap=function(EXP,N,Ckernel=0,Csubzones=0,HD=0){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); MAPS=NULL; 
if(Csubzones){Prob=numeric(N)}
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
MAPS[[i]]=createMAP(EXPv,1,Ckernel,MC=1)
if(Csubzones>0){Prob[i]=calcSubZones(MAPS[[i]],Csubzones)}}
if(Csubzones>0){P=createQuantiles(Prob); return(P)}
MAPS=seekMAPS(MAPS,N)
PlotMAP(MAPS[[1]],HD,main='5%ile') 
dev.new(); PlotMAP(MAPS[[3]],HD,main='95%ile')  
dev.new(); PlotMAP(MAPS[[2]],HD,main='mean') 
return(MAPS)}
MonteCarloMap(EXP[[1]],10^4)
Plot=function(H,main=NULL){
image(-H,zlim=c(-0.006,0),col=heat.colors(1000),main=main)}
PlotMAP=function(MAP,HD=0,main=NULL){
lev=c(0.1,0.25,0.5,1,1.5,2,3,4,5)
if(!HD){Plot(MAP,main); contour(MAP*1600,levels=lev,add=TRUE)}
if(HD){kernel_smoothMAP(MAP,11,lev)}}
MonteCarloMap(EXP[[1]],10^4)
MonteCarloMap(EXP[[1]],10^4,1)
MonteCarloMap(EXP[[1]],10^4,Csubzones=8)
createMAP=function(EXPv,flag=1,Ckernel=0,Csubzones=0,HD=0,MC=0){if(MC==0){EXPv=EXPv[2,]/100}
P=numeric(5)
P[1]=1-EXPv[1]
P[2]=EXPv[1]*EXPv[2]*EXPv[3]
P[3]=EXPv[1]*EXPv[2]*(1-EXPv[3])
P[4]=EXPv[1]*(1-EXPv[2])*EXPv[4]
P[5]=EXPv[1]*(1-EXPv[2])*(1-EXPv[4])
if(!flag){return(P*100)}
if(flag){
if(!Ckernel){MAP=maskTOT*P[1]+vents12*P[2]+vents3*P[3]+vfaults*P[4]+vfracts*P[5]}
if(Ckernel){MAP=maskTOT*P[1]+vkernel12*P[2]+vkernel3*P[3]+vfaults*P[4]+vfracts*P[5]}}
if(Csubzones>0){Prob=calcSubZones(MAP,Csubzones),return(Prob*100)}
if(flag==2){PlotMAP(MAP,HD)}
if(flag==1){return(MAP)}}
createMAP=function(EXPv,flag=1,Ckernel=0,Csubzones=0,HD=0,MC=0){if(MC==0){EXPv=EXPv[2,]/100}
P=numeric(5)
P[1]=1-EXPv[1]
P[2]=EXPv[1]*EXPv[2]*EXPv[3]
P[3]=EXPv[1]*EXPv[2]*(1-EXPv[3])
P[4]=EXPv[1]*(1-EXPv[2])*EXPv[4]
P[5]=EXPv[1]*(1-EXPv[2])*(1-EXPv[4])
if(!flag){return(P*100)}
if(flag){
if(!Ckernel){MAP=maskTOT*P[1]+vents12*P[2]+vents3*P[3]+vfaults*P[4]+vfracts*P[5]}
if(Ckernel){MAP=maskTOT*P[1]+vkernel12*P[2]+vkernel3*P[3]+vfaults*P[4]+vfracts*P[5]}}
if(Csubzones>0){Prob=calcSubZones(MAP,Csubzones);return(Prob*100)}
if(flag==2){PlotMAP(MAP,HD)}
if(flag==1){return(MAP)}}
createMAP(EXP[[1]],flag=2,Csubzones=8)
DM_ERF=function(EXP,Mexp,syntSEED){
pERF=ERFweights(Mexp,syntSEED)
DM=EXP[[1]]*0
for(i in 1:h){
DM=DM+EXP[[i]]*pERF[i]}
return(DM)}
DMerf=DM_ERF(EXP,Mexp,syntSEED)
plotEXPt(DMerf)
plotEXPt=function(EXP,i=0){if(i){EXP=EXP[[i]]}
plot(EXP[2,],seq(h,1),pch=16,xlim=c(0,100),ylim=c(0.5,h),col='red',main=c('Expert',i))
points(EXP[1,],seq(h,1),pch=1)
points(EXP[3,],seq(h,1),pch=1)
for(i in 1:h){
lines(seq(1,1000)/1000*100,rep(i,1000))
lines(seq(1,1000)/1000*(EXP[3,h+1-i]-EXP[1,h+1-i])+EXP[1,h+1-i],rep(i,1000),lwd=2)}}
plotEXPt(EXP,1)
plotEXPt(DMerf)
plotEXPt=function(EXP,i=0){if(i){EXP=EXP[[i]]}; if(!i){i='DM'}
plot(EXP[2,],seq(h,1),pch=16,xlim=c(0,100),ylim=c(0.5,h),col='red',main=c('Expert',i))
points(EXP[1,],seq(h,1),pch=1)
points(EXP[3,],seq(h,1),pch=1)
for(i in 1:h){
lines(seq(1,1000)/1000*100,rep(i,1000))
lines(seq(1,1000)/1000*(EXP[3,h+1-i]-EXP[1,h+1-i])+EXP[1,h+1-i],rep(i,1000),lwd=2)}}
plotEXPt(DMerf)
ERFweights(MexpS,syntSEED)
DM_ERF=function(EXP,Mexp,syntSEED){
pERF=ERFweights(Mexp,syntSEED)/100
DM=EXP[[1]]*0
for(i in 1:h){
DM=DM+EXP[[i]]*pERF[i]}
return(DM)}
plotEXPt(DMerf)
DMerf=DM_ERF(EXP,Mexp,syntSEED)
plotEXPt(DMerf)
syntEXP=function(Err,Inc,flag=1,i=1){
Bias=runif(k,-Err,Err)/100
Mexp=matrix(0,k,3)
Mexp[,2]=syntSEED*(1+Bias)
Mexp[,1]=Mexp[,2]*(1-Inc/100)
Mexp[,3]=Mexp[,2]*(1+Inc/100)
if(flag){plot(syntSEED,seq(1,k),pch=4,xlim=c(0,max(Mexp)),main=c('Expert',i),sub='Seed Questions' )
points(Mexp[,2],seq(1,10),col='red',pch=16)
points(Mexp[,1],seq(1,10),pch=1)
points(Mexp[,3],seq(1,10),pch=1)
for(i in 1:k){
lines(seq(1,1000)/1000*Mexp[,3],rep(i,1000))
lines(seq(1,1000)/1000*(Mexp[i,3]-Mexp[i,1])+Mexp[i,1],rep(i,1000),lwd=2)}}
return(round(Mexp))}
MexpS=createEXP(vErr,vInc,Q)
plotEXPt=function(EXP,i=0){if(i){EXP=EXP[[i]]}; if(!i){i='DM'}
plot(EXP[2,],seq(h,1),pch=16,xlim=c(0,100),ylim=c(0.5,h),col='red',main=c('Expert',i),sub='Target Questions')
points(EXP[1,],seq(h,1),pch=1)
points(EXP[3,],seq(h,1),pch=1)
for(i in 1:h){
lines(seq(1,1000)/1000*100,rep(i,1000))
lines(seq(1,1000)/1000*(EXP[3,h+1-i]-EXP[1,h+1-i])+EXP[1,h+1-i],rep(i,1000),lwd=2)}}
plotEXPt(DMerf)
?plot
plotEXPt=function(EXP,i=0){if(i){EXP=EXP[[i]]}; if(!i){i='DM'}
plot(EXP[2,],seq(h,1),pch=16,xlim=c(0,100),ylim=c(0.5,h),col='red',main=c('Expert',i),xlab='Target Questions',ylab=NULL)
points(EXP[1,],seq(h,1),pch=1)
points(EXP[3,],seq(h,1),pch=1)
for(i in 1:h){
lines(seq(1,1000)/1000*100,rep(i,1000))
lines(seq(1,1000)/1000*(EXP[3,h+1-i]-EXP[1,h+1-i])+EXP[1,h+1-i],rep(i,1000),lwd=2)}}
plotEXPt(DMerf)
plotEXPt=function(EXP,i=0){if(i){EXP=EXP[[i]]}; if(!i){i='DM'}
plot(EXP[2,],seq(h,1),pch=16,xlim=c(0,100),ylim=c(0.5,h),col='red',main=c('Expert',i),xlab='Responses',ylab='Target Questions')
points(EXP[1,],seq(h,1),pch=1)
points(EXP[3,],seq(h,1),pch=1)
for(i in 1:h){
lines(seq(1,1000)/1000*100,rep(i,1000))
lines(seq(1,1000)/1000*(EXP[3,h+1-i]-EXP[1,h+1-i])+EXP[1,h+1-i],rep(i,1000),lwd=2)}}
plotEXPt(DMerf)
plotEXPt=function(EXP,i=0){if(i){EXP=EXP[[i]]}; if(!i){i='DM'}
plot(EXP[2,],seq(h,1),pch=16,xlim=c(0,100),ylim=c(0.5,h),col='red',main=c('Expert',i),xlab='Target Responses',ylab='Target Questions')
points(EXP[1,],seq(h,1),pch=1)
points(EXP[3,],seq(h,1),pch=1)
for(i in 1:h){
lines(seq(1,1000)/1000*100,rep(i,1000))
lines(seq(1,1000)/1000*(EXP[3,h+1-i]-EXP[1,h+1-i])+EXP[1,h+1-i],rep(i,1000),lwd=2)}}
syntEXP=function(Err,Inc,flag=1,i=1){
Bias=runif(k,-Err,Err)/100
Mexp=matrix(0,k,3)
Mexp[,2]=syntSEED*(1+Bias)
Mexp[,1]=Mexp[,2]*(1-Inc/100)
Mexp[,3]=Mexp[,2]*(1+Inc/100)
if(flag){plot(syntSEED,seq(1,k),pch=4,xlim=c(0,max(Mexp)),main=c('Expert',i),xlab='Seed Responses', ylab='Seed Questions' )
points(Mexp[,2],seq(1,10),col='red',pch=16)
points(Mexp[,1],seq(1,10),pch=1)
points(Mexp[,3],seq(1,10),pch=1)
for(i in 1:k){
lines(seq(1,1000)/1000*Mexp[,3],rep(i,1000))
lines(seq(1,1000)/1000*(Mexp[i,3]-Mexp[i,1])+Mexp[i,1],rep(i,1000),lwd=2)}}
return(round(Mexp))}
MexpS=createEXP(vErr,vInc,Q)
plotEXPt(DMerf)
max_entropy=function(a,b,c,rA,rB){
x=runif(1)
y=runif(1,b,c)
if(x>0.95){y=runif(1,c,rB)}
if(x<0.5){y=runif(1,a,b)
if(x<0.05){y=runif(1,rA,a)}}
return(y)}
DM_linear=function(EXP,P=rep(1/k,k),N=10^5){
quan=matrix(0,h,3)
vexp=matrix(0,k,3)
for(i in 1:h){
for(j in 1:k){
vexp[j,]=EXP[[k]][,i]}
rA=min(vexp[,1])
rB=max(vexp[,3])
R=rB-rA
rA=rA-R/10
rB=rB+R/10
test=numeric(N)
for(j in 1:N){
estr=ceiling(runif(1)*k)
test[j]=max_entropy(vexp[estr,1],vexp[estr,2],vexp[estr,3],rA,rB)}
quan[i,]=createQuantiles(test,1)
print(i); flush.console()}
return(quan)}
DM_linear(EXP)
DM_linear=function(EXP,P=rep(1/k,k),N=10^5){
quan=matrix(0,h,3)
vexp=matrix(0,k,3)
for(i in 1:k){
for(j in 1:h){
vexp[j,]=EXP[[k]][,i]}
rA=min(vexp[,1])
rB=max(vexp[,3])
R=rB-rA
rA=rA-R/10
rB=rB+R/10
test=numeric(N)
for(j in 1:N){
estr=ceiling(runif(1)*k)
test[j]=max_entropy(vexp[estr,1],vexp[estr,2],vexp[estr,3],rA,rB)}
quan[i,]=createQuantiles(test,1)
print(i); flush.console()}
return(quan)}
DM_linear(EXP)
k
Q
DM_linear=function(EXP,P=rep(1/Q,Q),N=10^5){
quan=matrix(0,h,3)
vexp=matrix(0,Q,3)
for(i in 1:h){
for(j in 1:Q){vexp[j,]=EXP[[j]][,i]}
rA=min(vexp[,1]); rB=max(vexp[,3])
R=rB-rA; rA=rA-R/10; rB=rB+R/10
test=numeric(N)
for(j in 1:N){estr=ceiling(runif(1)*Q); test[j]=max_entropy(vexp[estr,1],vexp[estr,2],vexp[estr,3],rA,rB)}
quan[i,]=createQuantiles(test,1)
print(i); flush.console()}
return(quan)}
DM_linear(EXP)
?matrix
DM_linear=function(EXP,P=rep(1/Q,Q),N=10^5){
coln=c('5%ile','mean','95%ile')
quan=matrix(0,h,3,dimnames=list(coln,NULL))
vexp=matrix(0,Q,3)
for(i in 1:h){
for(j in 1:Q){vexp[j,]=EXP[[j]][,i]}
rA=min(vexp[,1]); rB=max(vexp[,3])
R=rB-rA; rA=rA-R/10; rB=rB+R/10
test=numeric(N)
for(j in 1:N){estr=ceiling(runif(1)*Q); test[j]=max_entropy(vexp[estr,1],vexp[estr,2],vexp[estr,3],rA,rB)}
quan[i,]=createQuantiles(test,1)
print(i); flush.console()}
return(quan)}
DM_linear=function(EXP,P=rep(1/Q,Q),N=10^5){
}
DM_linear=function(EXP,P=rep(1/Q,Q),N=10^5){
coln=c('5%ile','mean','95%ile')
quan=matrix(0,h,3,dimnames=list(coln,NULL))
vexp=matrix(0,Q,3)
for(i in 1:h){
for(j in 1:Q){vexp[j,]=EXP[[j]][,i]}
rA=min(vexp[,1]); rB=max(vexp[,3])
R=rB-rA; rA=rA-R/10; rB=rB+R/10
test=numeric(N)
for(j in 1:N){estr=ceiling(runif(1)*Q); test[j]=max_entropy(vexp[estr,1],vexp[estr,2],vexp[estr,3],rA,rB)}
quan[i,]=createQuantiles(test,1)
print(i); flush.console()}
return(quan)}
DM_linear(EXP)
DM_linear=function(EXP,P=rep(1/Q,Q),N=10^5){
coln=c('5%ile','mean','95%ile')
quan=matrix(0,h,3,dimnames=list(NULL,coln))
vexp=matrix(0,Q,3)
for(i in 1:h){
for(j in 1:Q){vexp[j,]=EXP[[j]][,i]}
rA=min(vexp[,1]); rB=max(vexp[,3])
R=rB-rA; rA=rA-R/10; rB=rB+R/10
test=numeric(N)
for(j in 1:N){estr=ceiling(runif(1)*Q); test[j]=max_entropy(vexp[estr,1],vexp[estr,2],vexp[estr,3],rA,rB)}
quan[i,]=createQuantiles(test,1)
print(i); flush.console()}
return(quan)}
DM_linear(EXP)
DM_linear(EXP,c(0.5,0.5,0,0))
syntSEED
syntSEED_fix=syntSEED
EXP
MexpS
MexpS_fix=MexpS
PlotMAP(vkernel12)
PlotMAP(vkernel3)
PlotMAP(vkernel12, HD=1)
PlotMAP(vkernel3); PlotMAP(vkernel3, HD=1)
PlotMAP(vfaults,HD=1)
PlotMAP(vfracts,HD=1)
syntSEED=round(runif(k,0,250))
syntSEED
MexpS=createEXP(vErr,vInc,Q)
ERFweights(MexpS,syntSEED)
MexpS
ERFweights
ERFsingle
ERFweights(MexpS,syntSEED)
ERFsingle(1,1,1,1)
syntSEED=round(runif(k,10,250)); syntSEED
ERFweights(MexpS,syntSEED)
MexpS=createEXP(vErr,vInc,Q)
ERFweights(MexpS,syntSEED)
plotEXPt(EXP,1)
createMAP(EXP[[1]],flag=0)
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
coln=c('EpoI_II','EpoIII','Faults','Fracts','Homog'); rown=c('5%ile','mean','95%ile')
a=matrix(0,3,5,dimnames=list(rown,coln))
for(i in 1:5){a[,i]=createQuantiles(coeff[,i])}
return(a)}
createMAP(EXP[[1]],flag=0)
createMAP=function(EXPv,flag=1,Ckernel=0,Csubzones=0,HD=0,MC=0){if(MC==0){EXPv=EXPv[2,]/100}
P=numeric(5)
P[1]=1-EXPv[1]
P[2]=EXPv[1]*EXPv[2]*EXPv[3]
P[3]=EXPv[1]*EXPv[2]*(1-EXPv[3])
P[4]=EXPv[1]*(1-EXPv[2])*EXPv[4]
P[5]=EXPv[1]*(1-EXPv[2])*(1-EXPv[4])
if(!flag){return(c(P[2:5],P[1])*100)}
if(flag){
if(!Ckernel){MAP=maskTOT*P[1]+vents12*P[2]+vents3*P[3]+vfaults*P[4]+vfracts*P[5]}
if(Ckernel){MAP=maskTOT*P[1]+vkernel12*P[2]+vkernel3*P[3]+vfaults*P[4]+vfracts*P[5]}}
if(Csubzones>0){Prob=calcSubZones(MAP,Csubzones);return(Prob*100)}
if(flag==2){PlotMAP(MAP,HD)}
if(flag==1){return(MAP)}}
createMAP(EXP[[1]],flag=0)
createMAP(EXP[[1]],flag=2)
createMAP(EXP[[1]],flag=2,Csubzones=8)
createMAP(EXP[[1]],flag=2,HD=1)
createMAP(EXP[[1]],flag=2,Ckernel=1)
MonteCarloCoeff(EXP[[1]],10^4)
MonteCarloMap(EXP[[1]],10^4)
outp=MonteCarloMap(EXP[[1]],10^4,1)
MonteCarloMap(EXP[[1]],10^4,Csubzones=8)
plotEXPt(DMerf)
DM_linear(EXP)
DM_linear(EXP,c(0.5,0.5,0,0))
DMerf
EXP[[1]]
MonteCarloCoeff(DMerf,10^4)
outp=MonteCarloMap(DMerf,10^4,1)
MonteCarloMap(DMerf,10^4,Csubzones=8)
outp=MonteCarloMap(DMerf,10^4,HD=1)
outp=MonteCarloMap(DMerf,10^4,HD=1)
outp=MonteCarloMap(DMerf,10^4)
PlotMAP=function(MAP,HD=0,main=NULL){
lev=c(0.1,0.25,0.5,1,1.5,2,3,4,5)
if(!HD){Plot(MAP,main); contour(MAP*1600,levels=lev,add=TRUE)}
if(HD){kernel_smoothMAP(MAP,11,lev,main)}}
kernel_smoothMAP=function(M,R,lev,main){
M=M/sum(M)
pencil=crea_tondi_SMAll(R)
n=ncol(M)
M2=matrix(0,5*n,5*n)
for(i in 1:(5*n)){for(j in 1:(5*n)){M2[i,j]=M[(i-1)%/%5+1,(j-1)%/%5+1]}}
M3=M2*0
for(i in R:(5*n-R+1)){for(j in R:(5*n-R+1)){
M3[i,j]=sum(M2[(i-R+1):(i+R-1),(j-R+1):(j+R-1)]*pencil)}}
Plot(M3,main=main); mM3=max(M3)
contour(M3*1600,add=T,levels=lev)}
lotMAP(vkernel12, HD=1)
PlotMAP(vkernel12, HD=1)
outp=MonteCarloMap(DMerf,10^4,HD=1)
q()