############################# # Programa_artigo: Simula??o# ############################# rm(list=ls(all=TRUE)) require(geoR) require(RandomFields) #Funções Utilizadas: source('/home/ajr/Área de Trabalho/TUDO/ESALQ/1MESTRADO/R/Dissertação/Simulação/funcaoana/lower2matrix.R') source('/home/ajr/Área de Trabalho/TUDO/ESALQ/1MESTRADO/R/Dissertação/Simulação/funcaoana/montaSigma.R') source('/home/ajr/Área de Trabalho/TUDO/ESALQ/1MESTRADO/R/Dissertação/Simulação/funcaoana/logpars2pars.R') source('/home/ajr/Área de Trabalho/TUDO/ESALQ/1MESTRADO/R/Dissertação/Simulação/funcaoana/geral.R') source('/home/ajr/Área de Trabalho/TUDO/ESALQ/1MESTRADO/R/Dissertação/Simulação/funcaoana/param.R') source('/home/ajr/Área de Trabalho/TUDO/ESALQ/1MESTRADO/R/Dissertação/Simulação/funcaoana/pred.R') LP <- list(nu1 = 0.9 , nu2 = 0.6, a = 1/7 , rho = 0.7 , sig1 = sqrt(0.7), sig2 = sqrt(0.5) , tau1 = sqrt(0.2), tau2 = sqrt(0.1) ) set.seed(7) sim<-1 nsample = 40 ntotal = 242 nnew = ntotal - 2*nsample gerado = matrix(NA,ncol=sim,nrow=ntotal ) variavel = matrix(NA,ncol=sim,nrow=nnew ) predi = matrix(NA,ncol=sim,nrow=2*nsample) ajustesim = matrix(NaN,sim,10) predicao = matrix(NA,ncol=sim,nrow=2*nsample ) krige.var = matrix(NA,ncol=sim,nrow=2*nsample ) llik = matrix(NA,ncol=sim,nrow=1) for (nsim in 1:sim) { gs <- cbind(runif(121, 0, 30), runif(121, 0, 30)) #distancia dM <- dist(as.matrix(gs), upper=T, diag=T) dV <- as.vector(dist(gs)) system.time(SIGMA <- montaSigma(LP, distPares=dV)) dM <- as.matrix(dist(gs, diag=T,up=T)); dimnames(dM) <- list(NULL, NULL) cholSIGMA <- chol(SIGMA) sim <- drop(crossprod(cholSIGMA, rnorm(dim(cholSIGMA)[1]))) rm(SIGMA) sim <- sim + c(rep(2, 121), rep(10, 121)) sim1 <- sim[1:(length(sim)/2)] sim2 <- sim[((length(sim)/2)+1):length(sim)] conj <- cbind(gs,sim1,sim2) amostra <- sample(seq(1:nrow(conj)),nsample,rep=FALSE) ###PREDIÇÃO conjunto.pred <- conj[amostra,] pred1 <- conjunto.pred[,3] pred2 <- conjunto.pred[,4] #UTILIZADO conjunto<-conj[-amostra,] var1<-conjunto[,3] var2<-conjunto[,4] lista<-list(n=length(c(var1,var2)), n1=length(var1), n2=length(var2), X=c(var1,var2), d=as.vector(dist(conjunto[,1:2]))) ini <-c(tau1=log(sqrt(0.4)), tau2=log(sqrt(0.3)), rho=log((0.75+1)/(1-(0.75))), nu1=1.3, nu2=0.5, sig2=log(sqrt(0.6)), a=log(500)) fixos <- c(mu1=5,mu2=1.2,sig1=log(sqrt(0.65))) ajuste <- optim(ini, geral, pars2=fixos, dados=lista,#method="BFGS", control=list(fnscale=-1, maxit=3000), conc=T, printpars =T,logp=T) ajuste$par <- logpars2pars(ajuste$par) pars.est<-list(a<-ajuste[[1]][[7]], rho<- (ajuste[[1]][[3]]), nu1<-ajuste[[1]][[4]], nu2<-ajuste[[1]][[5]], tau1<-(ajuste[[1]][[1]]), tau2<-(ajuste[[1]][[2]]), sig2<-(ajuste[[1]][[6]]), sig1<-1 ) names(pars.est) <- c('a','rho','nu1','nu2','tau1','tau2','sig2','sig1') dados<-list(n1=length(var1), n2=length(var2), X=c(var1,var2), d=as.vector(dist(conjunto[,1:2]))) parametros <- param(pars.est,dados) #ajustesim[nsim,1:10] <- parametros gerado[,nsim]<-c(sim1,sim2) variavel[,nsim]<-c(var1,var2) predi[,nsim]<-c(pred1,pred2) #} parametros<-list(mu1<-parametros[[1]], mu2<-parametros[[2]], sig1<-parametros[[3]], sig2<-parametros[[4]], tau1<-parametros[[5]], tau2<-parametros[[6]], nu1<-parametros[[7]], nu2<-parametros[[8]], a<-parametros[[9]], rho<-parametros[[10]] ) names(parametros) = c('mu1','mu2','sig1','sig2','tau1','tau2','nu1','nu2','a','rho') parametros ajustesim[nsim,1:10] <- as.numeric(parametros) ################## k = NULL k = as.data.frame(list(mu1,mu2,sig1,sig2,tau1,tau2,nu1,nu2,a,rho)) names(k) = c('mu1','mu2','sig1','sig2','tau1','tau2','nu1','nu2','a','rho') d = as.vector(dist(conjunto[,1:2])) S <- montaSigma(k, d) Sinv<-solve(S) object<-as.list(c(n1=length(conjunto[,3]), n2=length(conjunto[,4]), mu1=parametros$mu1, mu2=parametros$mu2, sigma1=parametros$sig1, sigma2=parametros$sig2, phi=parametros$a, kappa1=parametros$nu1, kappa2=parametros$nu2, rho=parametros$rho, tau1=parametros$tau1, tau2=parametros$tau2 )) geodata1<-as.geodata(conjunto,coords.col=1:2,data.col=3) geodata2<-as.geodata(conjunto,coords.col=1:2,data.col=4) ############# locations<-conjunto.pred[,1:2] predicao1<-pred(object,locations,variable.to.predict=1) #plot(conjunto.pred[,3],predição1$pred) pred.var1<-predicao1$pred krige.var1<-predicao1$krige.var predicao2<-pred(object,locations,variable.to.predict=2) #plot(conjunto.pred[,4],predição2$pred) pred.var2<-predicao2$pred krige.var2<-predicao2$krige.var predicao[,nsim]<-c(pred.var1,pred.var2) krige.var[,nsim]<-c(krige.var1,krige.var2) ##Loglik do modelo loglik <- geral(pars1=c(tau1=sqrt(tau1), tau2=sqrt(tau2), rho=rho, nu1=nu1, nu2=nu2, sig2=sqrt(sig2), a=a), pars2=c(mu1=mu1, mu2=mu2, sig1=sqrt(sig1)), dados=lista, conc=T,logpars=T, printpars=F) llik[,nsim]<-loglik } gerado variavel predi ajustesim predicao krige.var llik # # #Predi??o em um grid: # locations<- expand.grid(seq(min(conjunto[,1]),max(conjunto[,1]), by=0.2), # seq(min(conjunto[,2]),max(conjunto[,2]), by=0.2)) # # # # ##Erro Quadr?tico M?dio # EQM_v1=(predicao1$pred - var1.pred)^2 # EQM_v1 # summary(EQM_v1) # # EQM_v2=(predicao2$pred - var2.pred)^2 # EQM_v2 # summary(EQM_v2) # # ##Erro Absoluto # EA_v1=predicao1$pred - var1.pred # EA_v1 # summary(EA_v1) # # EA_v2=predicao2$pred - var2.pred # EA_v2 # summary(EA_v2) # # # #AIC # p= 10#n?mero de parametros # AIC<-2*(10-loglik) # AIC # # ##Intervalo de Predi??o # # Das vari?veis n?o utilizadas # # krige.var<-c(predicao1$krige.var,predicao2$krige.var) # predict<-c(predicao1$predict,predicao2$predict) # # IPinf<-predict-1.96*sqrt(krige.var) # IPinf # IPsup<-predict+1.96*sqrt(krige.var) # IPsup # # summary(IPinf,IPsup) # # variav<-c(var1.pred,var2.pred) # IO<-(variav>IPinf & variav