library(rgl) library(FNN) library(utils) library(markovchain) library(foreach) library(doParallel) library(Matrix) library(MASS) library(tikzDevice) library(RColorBrewer) my_palette <- colorRampPalette(c("blue", "green", "yellow", "red", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black", "black"))(n = 700) locCores<-24 horizon<-200 ###########CREATE ESTIMATION ERROR################# Atrue<-t(matrix(c(0.8,0.1,0.05, 0.15,0.8,0.15,0.05,0.1,0.8), nrow=3))#*0.25+diag(c(1,1,1))*0.75) clinktrue<-c(3,4,5)/2 clink2true<-c(5,4,3)/2 mcB<-new("markovchain", transitionMatrix=t(Atrue)) Xsim<-markovchainSequence(n=100,markovchain=mcB, t0="1") obsTrans<-createSequenceMatrix(Xsim) rowS<-outer(rowSums(obsTrans), c(1,1,1)) A<-t(obsTrans/rowS) a<-as.vector(t(A-diag(diag(A))-5*diag(c(1,1,1)))) a<-a[a!=-5] Amatrix<-function(a){ #Creates the A matrix from a vector #(forces columns to sum to 1, assumes entries in a are listed by row then column # without the diagonal) A<-matrix(data=0,nrow=3, ncol=3) n=1 for(i in 1:3){for(j in 1:3){ if(i!=j){A[i,j]=a[n] n=n+1} }} diag(A)<-1-colSums(A) A } #finds the approximate covariance matrix of a aErr<-matrix(data=0, nrow=6, ncol=6) state<-c(2,3,1,3,1,2) tCount<-as.vector(rowSums(obsTrans)[state]) diag(aErr)<-a*(1-a)/as.vector(rowSums(obsTrans)[state]) for(i in 1:6){for(j in 1:6){ if(i==j){ aErr[i,i]<-a[i]*(1-a[i])/tCount[i] }else if(state[i]==state[j]){ aErr[i,j]<--a[i]*a[j]/tCount[i] }}} aErr<-nearPD(aErr)$mat aPenMat<-as.matrix(solve(aErr)) Y<-rnorm(100,mean=clinktrue[as.integer(Xsim)]) Y2<-rnorm(100,mean=clink2true[as.integer(Xsim)]) clink<-c(aggregate(Y, list(Xsim), mean)$x, aggregate(Y2, list(Xsim), mean)$x) clinkErr<-1/(as.vector(summary(factor(Xsim)))) clinkErr<-diag(c(clinkErr, clinkErr)) cPenMat<-as.matrix(solve(clinkErr)) rm(state, tCount, rowS, obsTrans) ############SIMULATE DATA########################## X<-matrix(ncol=horizon, nrow=3) X[,1]<-c(1,0,0) Y<-0 Y2<-0 for(t in 2:horizon){ X[,t]<-rmultinom(1,1, Atrue%*%X[,t-1]) Y[t]<-rnorm(1, sum(clink[1:3]*X[,t])) Y2[t]<-rnorm(1, sum(clink[4:6]*X[,t])) } plot(Y,,col=c(1,2,3)%*%X) rescale<-function(x){x/sum(x)} #Basic filter Xhat<-matrix(ncol=horizon, nrow=3) Xhat[,1]<-c(1,1,1)/3 for(t in 2:horizon){ Xhat[,t]<-diag(dnorm(Y2[t], clink[4:6])*dnorm(Y[t], clink[1:3]))%*%A%*%Xhat[,t-1] Xhat[,t]<-rescale(Xhat[,t]) } XhatTrue<-matrix(ncol=horizon, nrow=3) XhatTrue[,1]<-c(1,1,1)/3 for(t in 2:horizon){ XhatTrue[,t]<-diag(dnorm(Y2[t], clink2true)*dnorm(Y[t], clinktrue))%*%Atrue%*%Xhat[,t-1] XhatTrue[,t]<-rescale(XhatTrue[,t]) } #par(mfrow=c(1,3)) #for(i in 1:3){ # plot(Xhat[i,],,'l') # lines(X[i,],,col='red') #} #par(mfrow=c(1,1)) logitmap<-function(x){ #returns logistic representation of components, fixing proportionality of first component to one #fails when components are zero log(x[-length(x)]/x[length(x)]) } alogit<-function(x){ y<-c(exp(x),1) y/sum(y) } #plot on simplex simplex.y <- function(x) { return( sqrt(0.75) * x[3] / sum(x) ) } simplex.x <- function(x) { return( (x[2] + 0.5 * x[3]) / sum(x) ) } bounds<-matrix(c(1,0,0,0,1,0,0,0,1,1,0,0), ncol=4) bounds<-list(x=apply(bounds, 2,simplex.x), y=apply(bounds, 2,simplex.y)) # pdf("Filterplot.pdf", width=14, height=6) textwidth<-4.7 tikz('Filterplot.tex', standAlone=F, width=textwidth, height=0.45*textwidth) par(mfrow=c(1,2), mar=c(1.3,0,1,1)) #image(plotx,ploty, log(1e-5+plotz), col=my_palette) XhatTrueplot<-list(x=apply(XhatTrue[,1:100], 2,simplex.x), y=apply(XhatTrue[,1:100], 2,simplex.y)) plot(bounds,, 'l', col='red', xlab="", ylab="", axes=F, xlim=c(-0.13,1.13), ylim=c(-0.08, 0.95)) title(sub="True Parameters", line=0) lines(XhatTrueplot, col='grey') points(XhatTrueplot, cex=0.8) text(-0.08,0.0, "$e_1$", cex=1) text(1.08,0.0, "$e_2$", cex=1) text(0.5,0.95, "$e_3$", cex=1) Xhatplot<-list(x=apply(Xhat[,1:100], 2,simplex.x), y=apply(Xhat[,1:100], 2,simplex.y)) plot(bounds,, 'l', col='red', xlab="", ylab="", axes=F, xlim=c(-0.13,1.13), ylim=c(-0.08, 0.95)) title(sub="Statistical Parameters", line=0) lines(Xhatplot, col='grey') points(Xhatplot, cex=0.8) text(-0.08,0.0, "$e_1$", cex=1) text(1.08,0.0, "$e_2$", cex=1) text(0.5,0.95, "$e_3$", cex=1) dev.off() par(mfrow=c(1,1)) lXhat<-apply(Xhat, 2, logitmap) plot(t(apply(Xhat, 2, logitmap))) # pdf("Filterplot.pdf", width=14, height=6) pptz<-t(apply(Xhat, 2, logitmap)) textwidth<-4.7 tikz('LogitFilter.tex', standAlone=F, width=textwidth, height=0.45*textwidth) par(mfrow=c(1,1), mar=c(5,5,1,1)) plot(pptz,, 'l', col='grey', xlab="$\\log(p_1/p_3)$", ylab="$\\log(p_2/p_3)$") title(sub="Statistical Parameters", line=-1) points(pptz) dev.off() ################# Derivatives of filter steps######### #Dl/Dp dlambda<-function(p){ D<-matrix(data=0,nrow=2, ncol=3) D[,3]<- -1/p[3] D[,-3]<- diag(1/p[-3]) D } #Drescale/Dp dR<-function(p){ diag(c(1,1,1)*1/sum(p))-outer(p, c(1,1,1)/sum(p)^2) } #DlambdaInv/Dp (as function of l) dlambdainv<-function(l){ p<-c(exp(l),1) dR(p)%*%(diag(p)[,-3]) } #Dp/Da (does not include correction) dApda<-function(pm){ #looks weird but gives right structure t(matrix(data=c( pm[2], pm[3], -pm[1], 0, -pm[1], 0, -pm[2], 0, pm[1], pm[3], 0, -pm[2], 0, -pm[3], 0, -pm[3], pm[1], pm[2]) , nrow=6, ncol=3)) } #Dp/Dclink (does not include prediction, does include rescaling) dRCpdc<-function(p,clink,y){ phat<-dnorm(y[2], clink[4:6])*dnorm(y[1], clink[1:3])*p D1<-outer(phat,c(1,1,1,1,1,1)) pd<-c(y[1]*c(1,1,1), y[2]*c(1,1,1))-clink D2<- t(matrix(data=c( pd[1],0,0,pd[4],0,0, 0,pd[2],0,0,pd[5],0, 0,0,pd[3],0,0,pd[6] ) , nrow=6, ncol=3)) dR(phat)%*%(D2*D1) } #save.image("FNNfittedStatW.RData") ### Penalty calculation####### D2kappa<-list() D2kappa[[1]]<-diag(c(1,1)) abest<-list() cbest<-list() for(t in 2:horizon){ #prediction step d2k<-D2kappa[[t-1]] C<-diag(dnorm(Y2[t], clink[4:6])*dnorm(Y[t], clink[1:3])) DRescalLambda<-dlambda(rescale(C%*%A%*%alogit(lXhat[,t-1])))%*% dR(as.vector(C%*%A%*%alogit(lXhat[,t-1]))) psir<-DRescalLambda%*%C%*%A%*%dlambdainv(lXhat[,t-1]) psiA<-DRescalLambda%*%C%*%dApda(Xhat[,t-1]) psiC<-DRescalLambda%*%C%*%dRCpdc(as.vector(A%*%Xhat[,t-1]),clink,c(Y[t],Y2[t])) psiAC<-cbind(psiA, psiC) psirI<-solve(psir) psirIpsiAC<-psirI%*%psiAC Sigma<-solve(bdiag(aPenMat, cPenMat) + t(psirIpsiAC)%*%d2k%*%psirIpsiAC) B<-t(psirI)%*%d2k%*%psirIpsiAC H<-t(psirI)%*%d2k%*%psirI D2kappa[[t]]<-as.matrix(nearPD(H-B%*%Sigma%*%t(B))$mat) #Could also use simplified equation for inverse, but this wouldn't give approximate starting points #for regression methods as easily. #Needs check for positive definiteness in case of numerical error QapOptimizer<-Sigma%*%t(B) abest[[t]]<-as.matrix(QapOptimizer[1:6,]) cbest[[t]]<-as.matrix(QapOptimizer[7:12,]) } kappaQ<-function(lp,k){ as.numeric(t(lp-lXhat[,k])%*%D2kappa[[k]]%*%(lp-lXhat[,k])) } ###########SET UP PENALTIES################################################ #clinkAlt<-clink+c(1,1,1) #clink2Alt<-clink2+c(1,1,1) #AAlt<-A*0.9+ diag(c(1,1,1))*0.1 neighb=25 #number of neighbours to use searchpar<-3000 #maximum number of approximation points optimsteps<-7000 #maximum number of optimization steps bigval<-1e5 #create basic set of reference points in the simplex logbase<-matrix(2*rnorm(2*searchpar), nrow=2) penalty<-function(aAlt, clinkAlt){ t(aAlt-a)%*%aPenMat%*%(aAlt-a) + t(clinkAlt-clink)%*%cPenMat%*%(clinkAlt-clink) } #create penalty function wrapper kappa<-function(lp,k){ out<-as.numeric(knn.reg(kappafit[[k]]$pred, lp, kappafit[[k]]$kappa, k=neighb)$pred-error[k]) if(is.na(out)){1e5}else{out} } #initialize penalty kappafit<-list() kappafit[[1]]<-list(pred=t(cbind(logbase, c(0,0))), kappa=c(logbase[1,]^2+logbase[2,]^2,0)) error<-knn.reg(kappafit[[1]]$pred, lXhat[,1], kappafit[[1]]$kappa, k=neighb)$pred runmax<-1 cl<-makeCluster(locCores) registerDoParallel(cl) timer<-Sys.time() for(t in (runmax+1):horizon){ #find perturbation points around which to scan Xhatpert<-logbase+outer(lXhat[,t], rep.int(1,searchpar)) kappacalc<-0 kappacalc<-foreach(n=(1:dim(Xhatpert)[2]), .combine='c') %dopar% { library('FNN') kappac<-function(par){ AAlt=Amatrix(par[1:6]) clinkAlt<-par[7:12] p<-par[13:14] optpen<-sum(abs(alogit(Xhatpert[,n])-rescale(diag(dnorm(Y2[t], clinkAlt[4:6])*dnorm(Y[t], clinkAlt[1:3]))%*%AAlt%*%alogit(p))))+ sum(pmax(-AAlt,0)) optpen<-(optpen)*bigval kappa(p, t-1)+penalty(par[1:6], par[7:12])+optpen } astart<-abest[[t]]%*%(Xhatpert[,n]-lXhat[,t])+a cstart<-cbest[[t]]%*%(Xhatpert[,n]-lXhat[,t])+clink pstart<-logitmap(pmax(solve(diag(dnorm(Y2[t], cstart[4:6])*dnorm(Y[t], cstart[1:3]))%*%Amatrix(astart), alogit(Xhatpert[,n])),1e-4)) optim(c(astart,cstart, pstart), kappac, control=list(maxit=optimsteps, reltol=1e-6))->opt final<-kappa(opt$par[13:14], t-1)+penalty(opt$par[1:6], opt$par[7:12]) if(abs(opt$value-final)/final>0.1){ bigval<-1e7 optim(opt$par, kappac, control=list(maxit=optimsteps, reltol=1e-6))->opt final<-kappa(opt$par[13:14], t-1)+penalty(opt$par[1:6], opt$par[7:12]) bigval<-1e5 } final } #plot3d(apply(Xhatpert, 2,simplex.x), apply(Xhatpert, 2,simplex.y), kappacalc, col='red', size=3) kappafit[[t]]<-list(pred=t(Xhatpert), kappa=kappacalc) error[t]<-knn.reg(kappafit[[t]]$pred, lXhat[,t], kappafit[[t]]$kappa, k=neighb)$pred timer[t]<-Sys.time() print(t) runmax<-max(runmax, t) #save.image("BothWUFSFly.RData") } stopCluster(cl) #plot3d(Xhatpert[1,], Xhatpert[2,], kappafit$fitted-error, col='red', size=3) #plot3d(Xhatpert[1,], Xhatpert[2,], kappafit[[t]]$kappa-error, col='red', size=3) #good<-(kappafit[[t]]$kappa<30) #plot3d(Xhatpert[1,good], Xhatpert[2,good], kappafit[[t]]$kappa[good]-error, col='red', size=3) #save.image("BothApproxWUFS.RData") ############PLOT#################### #plot the penalty on the simplex # plotx<-apply(basepoints, 2,simplex.x) # ploty<-apply(basepoints, 2,simplex.y) # plotz<-apply(basepoints, 2, kappa, k=t) # plot3d(plotx, ploty, plotz, col='red', size=3) #Xhatpertplot<-list(x=apply(Xhatpert, 2,simplex.x), y=apply(Xhatpert, 2,simplex.y)) #plot(bounds,, 'l', col='red') #points(Xhatpertplot) #pdf("Stat.pdf") # par(mfrow=c(round(sqrt(horizon)), round(horizon/round(sqrt(horizon))))) # par(mar=c(1,1,3,1)) # plottimes=round(sqrt(horizon))*round(horizon/round(sqrt(horizon))) # # for(t in 1:min(plottimes, horizon, runmax)){ pix<-100 plotz2<-plotz<-matrix(NA,nrow=pix, ncol=pix) plotx<-ploty<-(1:pix)/(1+pix) for(i in 1:pix){ for(j in 1:pix){ p3<-sqrt(4/3)*ploty[j] p2<-plotx[i]-p3/2 p1<-1-p2-p3 if (min(p1,p2,p3)<0){plotz[i,j]<-NA} else { plotz[i,j]<-kappaQ(logitmap(c(p1,p2,p3)),t) plotz2[i,j]<-kappa(logitmap(c(p1,p2,p3)),t) #plotz[i,j]<-sum(logitmap(c(p1,p2,p3))^2) #if(plotz[i,j]>1e3){plotz[i,j]<-NA} } } } # par(mfrow=c(1,2), mar=c(2,1,1,1)) # #image(plotx,ploty, log(1e-5+plotz), col=my_palette) # contour(plotx,ploty,plotz-min(plotz, na.rm=T), zlim=c(0, 10), # axes=F, xlab="", ylab="", levels=(1:20)/2) # title(main="Quadratic Approx.", line=-2) # title(sub=paste("Time t =", t), line=0) # points(simplex.x(Xhat[,t]),simplex.y(Xhat[,t]), cex=1, col='purple' ) # points(simplex.x(XhatTrue[,t]),simplex.y(XhatTrue[,t]), cex=2, col='purple', pch="*" ) # lines(bounds,,col='red') # # contour(plotx,ploty,plotz2-min(plotz2, na.rm=T), zlim=c(0, 10), # axes=F, xlab="", ylab="", levels=(1:20)/2) # title(main="KNN-Regression Approx.", line=-2) # title(sub=paste("Time t =", t), line=0) # points(simplex.x(Xhat[,t]),simplex.y(Xhat[,t]), cex=1, col='purple' ) # points(simplex.x(XhatTrue[,t]),simplex.y(XhatTrue[,t]), cex=2, col='purple', pch="*" ) # lines(bounds,,col='red') # par(mfrow=c(1,1), mar=c(5.1,4.1,4.1,2.1)) textwidth<-4.7 tikz(paste('penalty',t,".tex", sep=""), standAlone=F, width=textwidth, height=0.45*textwidth) par(mfrow=c(1,2), mar=c(1.3,0,1,1)) #image(plotx,ploty, log(1e-5+plotz), col=my_palette) contour(plotx,ploty,plotz-min(plotz, na.rm=T), zlim=c(0, 10), axes=F, xlab="", ylab="", levels=(1:10), labcex=0.4, xlim=c(-0.13,1.13), ylim=c(-0.08, 0.95)) title(sub=paste("Quadratic Approx., $t=", t, "$"), line=0) points(simplex.x(Xhat[,t]),simplex.y(Xhat[,t]), cex=1, col='purple' ) points(simplex.x(XhatTrue[,t]),simplex.y(XhatTrue[,t]), cex=0.8, col='purple', pch="*" ) lines(bounds,,col='red') text(-0.08,0.0, "$e_1$", cex=1) text(1.08,0.0, "$e_2$", cex=1) text(0.5,0.95, "$e_3$", cex=1) contour(plotx,ploty,plotz2-min(plotz2, na.rm=T), zlim=c(0, 10), axes=F, xlab="", ylab="", levels=(1:10), labcex=0.4, xlim=c(-0.13,1.13), ylim=c(-0.08, 0.95)) title(sub=paste("KNN-Regression Approx., $t=", t, "$"), line=0) points(simplex.x(Xhat[,t]),simplex.y(Xhat[,t]), cex=1, col='purple' ) points(simplex.x(XhatTrue[,t]),simplex.y(XhatTrue[,t]), cex=0.8, col='purple', pch="*" ) lines(bounds,,col='red') text(-0.08,0.0, "$e_1$", cex=1) text(1.08,0.0, "$e_2$", cex=1) text(0.5,0.95, "$e_3$", cex=1) dev.off() } # # #dev.off() # save.image("FNNfittedStat.RData") ############# Plot P(X_t=e_1) Under nonlinear expectations ####### cl<-makeCluster(locCores) registerDoParallel(cl) risk<-20 Px<-foreach(n=(1:min(plottimes, horizon)), .combine='c') %dopar% { library('FNN') penexp<-function(r){ p<-exp(r) -risk*p[1]/(1+sum(p))+kappaQ(r,n) } -optim(lXhat[,n], penexp,control=list(maxit=optimsteps))$value/risk } Pmx<-foreach(n=(1:min(plottimes, horizon)), .combine='c') %dopar% { library('FNN') penexp<-function(r){ p<-exp(r) risk*p[1]/(1+sum(p))+kappaQ(r,n) } optim(lXhat[,n], penexp,control=list(maxit=optimsteps))$value/risk } stopCluster(cl) Pmx10<-Pmx Px10<-Px Pmx5<-Pmx Px5<-Px plot(Xhat[1,1:100],,'l', xlim=c(0,100), ylim=c(0,1), ylab='P(Xt=e1)', xlab='Time') #lines(Px, col='red') #lines(Pmx, col='blue') polygon(c(1:horizon,rev(1:horizon)), c(Pmx, rev(Px)), col='skyblue') lines(Xhat[1,1:100], col='blue') textwidth<-4.7 tikz(paste("PX1.tex", sep=""), standAlone=F, width=textwidth, height=0.6*textwidth) plot(Xhat[1,1:100],,'l', xlim=c(1,20), ylim=c(0,1), ylab='$P(X_t=e_1)$', xlab='Time') mar=c(1.3,0,1,1) #lines(Px, col='red') #lines(Pmx, col='blue') polygon(c(1:horizon,rev(1:horizon)), c(Pmx, rev(Px)), col='lightgrey') polygon(c(1:horizon,rev(1:horizon)), c(Pmx10, rev(Px10)), col='grey') polygon(c(1:horizon,rev(1:horizon)), c(Pmx5, rev(Px5)), col='darkgrey') lines(Xhat[1,1:100]) points(XhatTrue[1,1:100], pch="*") dev.off() #######Entropy comparisons########## EntQ<-D2kappa EntMean<-0 EntSd<-0 EntMAD<-0 for(n in 1:horizon){ p<-Xhat[,n] #EntQ[[n]]<-diag(c(1/p[1]^3, 1/p[2]^3))+1/p[3]^3 EntQ[[n]]<-t(dlambdainv(logitmap(p)))%*%diag(1/p)%*%dlambdainv(logitmap(p)) RatMat<-D2kappa[[n]]/EntQ[[n]] EntMean[n]<-mean(RatMat) EntSd[n]<-sd(RatMat) EntMAD[n]<-mean(abs(RatMat-mean(RatMat))) } t<-70 EntT<-0 kapQ<-0 for(i in 1:3000){ p<-alogit(kappafit[[t]]$pred[i,]) EntT[i]<-sum(p*log(p/Xhat[,t])) kapQ[i]<-kappaQ(kappafit[[t]]$pred[i,],t) } textwidth<-4.7 tikz(paste("entropy.tex", sep=""), standAlone=F, width=textwidth, height=0.5*textwidth) par(mfrow=c(1,2), mar=c(4,5,1,1)) plot(EntT, kappafit[[t]]$kappa, xlab="$H(p|p^*_{70})$", ylab="$\\kappa(p,70)$", cex=0.6) plot(EntT, kappafit[[t]]$kappa, xlim=c(0,0.1), ylim=c(0,30), xlab="$H(p|p^*_{70})$", ylab="$\\kappa(p,70)$", cex=0.6) dev.off() plot(EntT, kapQ, xlim=c(0,0.1), ylim=c(0,30)) plot(kapQ,kappafit[[t]]$kappa, xlim=c(0,30), ylim=c(0,30)) ####Plot error of lXhat against minimum point### mXhat<-matrix(nrow=2, ncol=horizon) miner<-0 for(n in 2:horizon){ penexp<-function(r){ kappa(r,n) } test<-optim(lXhat[,n], penexp,control=list(maxit=optimsteps)) mXhat[,n]<-test$par miner[n]<--test$value } textwidth<-4.7 tikz(paste("kappa0.tex", sep=""), standAlone=F, width=textwidth, height=0.5*textwidth) par(mar=c(5,5,1,1)) plot(miner,, log='y', ylab="$\\kappa(p^*_t,t)$", xlab="Time") dev.off() ### Plot penalty evaluated at true parameter filter #### truePen<-0 for(t in 1:horizon){ truePen[t]<-kappaQ(logitmap(XhatTrue[,t]),t) } plot(truePen) #########CREATE ANIMATION###### library("animation") # # cl<-makeCluster(4) # registerDoParallel(cl) # # saveLatex({ # par(mar = c(3, 3, 1, 0.5), mgp = c(2, 0.5, 0), tcl = -0.3, cex.axis = 0.8, # cex.lab = 0.8, cex.main = 1) # for(t in 1:200){ # pix<-50 # plotz<-matrix(NA,nrow=pix, ncol=pix) # plotx<-ploty<-(0:pix)/pix # #plotz<-foreach(i=1:pix, .combine='rbind') %dopar% { # # library('FNN') # for(i in 1:pix){ # #pzcol<-0 # for(j in 1:pix){ # p3<-sqrt(4/3)*ploty[j] # p2<-plotx[i]-p3/2 # p1<-1-p2-p3 # if (min(p1,p2,p3)<0){plotz[i,j]<-Inf}#pzcol[j]<-Inf} # else { # plotz[i,j]<-kappa(logitmap(c(p1,p2,p3)),t) # # p<-kappa(logitmap(c(p1,p2,p3)),t) # #if(p>1e3){pzcol[j]<-Inf}else{pzcol[j]<-p} # } # } # #pzcol # } # image(plotx,ploty,plotz, zlim=c(0, 20), col=my_palette, main=t, axes=F, xlab="", ylab="") # points(simplex.x(Xhat[,t]),simplex.y(Xhat[,t]), cex=2, col='purple' ) # points(simplex.x(X[,t]*0.95+0.05*c(1,1,1)),simplex.y(X[,t]*0.95+0.05*c(1,1,1)), cex=3, pch='*', col='purple' ) # lines(bounds,,col='red') # # print(t) # } # }, img.name = "Anim3", ani.opts = "controls,loop,width=0.95\\textwidth", # latex.filename = ifelse(interactive(), "STAMPy2.tex", ""), # interval = 0.1, nmax = horizon, ani.dev = "pdf", ani.type = "pdf", ani.width = 7, # ani.height = 7, documentclass = paste("\\documentclass{article}", # "\\usepackage[papersize={7in,7in},margin=0.3in]{geometry}", # sep = "\n")) # stopCluster(cl) ## the PDF graphics output is often too large because it is ## uncompressed; try the option ani.options('pdftk') or ## ani.options('qpdf') to compress the PDF graphics; see ?pdftk or ## ?qpdf and ?ani.options #save.image("medCaseGood.RData") #try plotting on logistic space. pix<-100 plotz<-matrix(NA,nrow=pix, ncol=pix) cl<-makeCluster(locCores) registerDoParallel(cl) for(t in 1:min(horizon, plottimes, runmax)){ plotx<-5*(1:pix-pix/2)/(pix+1)+lXhat[1,t] ploty<-5*(1:pix-pix/2)/(pix+1)+lXhat[2,t] plotz<-foreach(i=1:(pix), .combine='rbind') %dopar% { library('FNN') pzcol<-0 for(j in 1:(pix)){ p<-kappa(c(plotx[i], ploty[j]),t) if(p>1e3){pzcol[j]<-Inf}else{pzcol[j]<-p} } as.vector(pzcol) } plotz2<-foreach(i=1:(pix), .combine='rbind') %dopar% { pzcol<-0 for(j in 1:(pix)){ p<-kappaQ(c(plotx[i], ploty[j]),t) if(p>1e3){pzcol[j]<-Inf}else{pzcol[j]<-p} } as.vector(pzcol) } #for(i in 1:(pix)){ # for(j in 1:pix){ # plotz[i,j]<-kappaQ(c(plotx[i], ploty[j]),t)-kappa(c(plotx[i], ploty[j]),t) # } #} par(mfrow=c(1,2)) contour(plotx,ploty,plotz-min(plotz),col=my_palette, zlim=c(0, 5), main=t) points(lXhat[1,t],lXhat[2,t], cex=2, col='purple', pch=20) contour(plotx,ploty,plotz2-min(plotz2),col=my_palette, zlim=c(0, 5),main=t) points(lXhat[1,t],lXhat[2,t], cex=2, col='purple', pch=20) par(mfrow=c(1,1)) } stopCluster(cl)