/**************************/ /* Load R-libraries used */ /**************************/ library("boot") library("class") library("foreign") library("MASS") library("sp") library("maptools") library("msm") library("nnet") library("lattice") library("Matrix") library("spatial") library("SparseM") library("spam") library("fields") library("tripack") library("tools") library("deldir") library("nlme") library("coda") library("spdep") /********************************************/ /* Simulate noisy image on a grid */ /********************************************/ /* define Moore's Basis functions: s=range of angular distance: q/2pi n=number of basis functions d=Dmax */ b<-function(s,n,d){ B<-(d/2)*(((sin(pi*n-2*pi*d*s)/(pi*n-2*pi*d*s))-(sin(pi*n+2*pi*d*s)/(pi*n+2*pi*d*s)))) return(B) } /*generate noiseless image intensities on a grid s=angular distance:q/2pi grid size: 100x100 beam-center: (50,50) s range: (0.001,0.075) or q range: (0.006,0.47) a<-c(5.281465e-05, 4.178242e-05, 1.686120e-05 ,5.454271e-06, 3.379092e-06, 3.007528e-06, 2.055125e-06, 8.427587e-07) # 8 Moore's coefficients d<-58.5 # Dmax */ x1<-seq(0.001,sqrt(2)*0.075,length=100) y1<-seq(0.001,sqrt(2)*0.075,length=100) # x1, y1 defines the grid corresponding to the s range S<-matrix(0,ncol=length(x1), nrow=length(y1)) # matrix to store the s distance from beam-center K<-matrix(0,ncol=length(x1), nrow=length(y1)) I_data<-matrix(0,ncol=length(x1), nrow=length(y1)) # Noiseless intensity matrix for (i in 1:length(x1)) { for(j in 1:length(y1)) { S[i,j]<-sqrt((x1[i]-round(x1[50],3))^2+(y1[j]-round(x1[50],3))^2) K[i,j]<-(a[1]*b(S[i,j],1,d)+a[2]*b(S[i,j],2,d)+a[3]*b(S[i,j],3,d)+a[4]*b(S[i,j],4,d)+a[5]*b(S[i,j],5,d)+a[6]*b(S[i,j],6,d)+a[7]*b(S[i,j],7,d)+a[8]*b(S[i,j],8,d)) I_data[i,j]<-((K[i,j])/S[i,j])*100000 } } /*generate the noisy image intensities */ nbd<-cell2nb(100,100,type="queen") # setup spatial neighorhood structure of image pixels: queen type nbd chosen N<-nb2listw(nbd,style="W") NN=similar.listw(N) W=as_dsTMatrix_listw(NN) A=as_dsCMatrix_IrW(W,rho=0.5) AA=as(A,"CsparseMatrix") Dinv=as(Diagonal(x=1/I_data),"CsparseMatrix") X<-(Dinv%*%AA) XpX<-t(X)%*%X error=rnorm(10000,mean=0,sd=1) Xpe<-t(X)%*%error sol<-as.vector(solve(XpX,Xpe)) E<-matrix(sol,nrow=100) #error matrix data=I_data+E # final noisy image /***************************************/ /* Analysis of a pixilated image */ /***************************************/ /* likelihood function for the spatial independent case */ lpostf<-function(a1,a2,a3,a4,a5,a6,a7,a8,y=data,cxx,cyy,d){ cx=floor(cxx) cy=floor(cyy) del_cx=cxx-cx del_cy=cyy-cy I<-rep(0,10000) res<-rep(0,10000) S<-rep(0,10000) Kfunc<-rep (0,10000) x1<-seq(0.001,sqrt(2)*0.075,length=100) y1<-seq(0.001,sqrt(2)*0.075,length=100) g=expand.grid(x=x1,y=y1) S<-sqrt((g\$x-x1[cx]+0.000001-del_cx*0.001061273)^2+(g\$y-y1[cy]+0.000001-del_cy*0.001061273)^2) a<-c(a1,a2,a3,a4,a5,a6,a7,a8)/100000 Kfunc<-(a[1]*b(S,1,d)+a[2]*b(S,2,d)+a[3]*b(S,3,d)+a[4]*b(S,4,d)+a[5]*b(S,5,d)+a[6]*b(S,6,d)+a[7]*b(S,7,d)+a[8]*b(S,8,d)) U<-ifelse((Kfunc > 0)==T,Kfunc,0.00001) I<-(U*100000)/S res<-(y-I) Vhat=as.vector(I) llike=(sum(((res^2)/(2*Vhat)),na.rm=T))+0.5*sum(log(Vhat),na.rm=T) prior_part<-sum(log(a))+sum((0.22*(log(a[1:8])-1)^2))+2*d-109*log(d)+(((cxx-50)^2)/8)+(((cyy-50)^2)/8) logpost<-(-(prior_part+llike)) return(logpost) } /* MCMC updates of all parameters except the autocorrelation NUMIT: number of iterations of MCMC s,s1,s2,s3: std errors used to increment parameter values for MCMC updates */ MHsampler = function(NUMIT,s,s1,s2,s3) { n = length(as.vector(data)) cat("n=",n,"\n") ## NUMIT matrix to store Markov chain values ## each row corresponds to one of 11 parameters in order: a's,cx,cy,d ## each column corresponds to a single state of the Markov chain mchain = matrix(NA,15,NUMIT) acc_a12d=0 # acceptance rates for parameter updates acc_a3=0 acc_a4=0 acc_a5=0 acc_a6=0 acc_a7=0 acc_a8=0 acc_cx=0 acc_cy=0 ## starting values for Markov chain log_curr<-lpostf(5.26,4.18,1.70,0.55,0.34,0.30,0.21,0.09,data,50.00,49.99,58.71) # initial call to evaluate logl rho1d=-0.99 ### pairwise correlations for joint updates of a1,a2 and Dmax rho2d= 0.93 rho12=-0.93 mchain[,1] =c(5.26,4.18,1.70,0.55,0.34,0.30,0.21,0.09,50.00,49.99,58.71,log_curr,rho1d,rho2d,rho12) #initial parameter values mu12d=c(mchain[1,1],mchain[2,1],mchain[11,1]) sa1=0.03 sa2=0.007 Sigma12d=matrix(c(sa1^2,rho12*sa1*sa2,rho1d*sa1*s2, rho12*sa1*sa2,sa2^2,rho2d*sa2*s2,rho1d*sa1*s2,rho2d*sa2*s2,s2^2),3,3) ### 3 D correlation matrix for joint updates for (i in 2:NUMIT) { ## most upto date state for each parameter curra1 = mchain[1,i-1] curra2 = mchain[2,i-1] curra3 = mchain[3,i-1] curra4 = mchain[4,i-1] curra5 = mchain[5,i-1] curra6 = mchain[6,i-1] curra7 = mchain[7,i-1] curra8 = mchain[8,i-1] currcx = mchain[9,i-1] currcy = mchain[10,i-1] currd = mchain[11,i-1] cat(" Iteration:",i) ## sample from full conditional distribution of theta ( joint Metropolis-Hastings update) propa12d = mvrnorm(1,mu12d,Sigma12d) # draw one sample at random ) propa1=propa12d[1] propa2=propa12d[2] propd=propa12d[3] ## Metropolis accept-reject step (in log scale) log_prop<-lpostf(propa1,propa2,curra3,curra4,curra5,curra6,curra7,curra8,data,currcx,currcy,propd) logMHratio1=log_prop-log_curr logalpha1 = min(0,logMHratio1) # alpha = min(1,MHratio) if ((log(runif(1))<(log(1)+logalpha1)) ) # accept if unif(0,1) 0){ rho1d=cor(mchain[1,1:batch[batch_index]],mchain[11,1:batch[batch_index]]) ### update the correltion matrix of the parameters a11, a12 and Dmax at regular interval rho2d=cor(mchain[2,1:batch[batch_index]],mchain[11,1:batch[batch_index]]) rho12=cor(mchain[1,1:batch[batch_index]],mchain[2,1:batch[batch_index]]) sa1=sqrt(var(mchain[1,1:batch[batch_index]])) sa2=sqrt(var(mchain[2,1:batch[batch_index]])) s2=sqrt(var(mchain[11,1:batch[batch_index]])) mua1=mean(mchain[1,1:batch[batch_index]]) mua2=mean(mchain[2,1:batch[batch_index]]) mud=mean(mchain[11,1:batch[batch_index]]) mu12d=c(mua1,mua2,mud) Sigma12d=matrix(c(sa1^2,rho12*sa1*sa2,rho1d*sa1*s2,rho12*sa1*sa2,sa2^2,rho2d*sa2*s2,rho1d*sa1*s2,rho2d*sa2*s2,s2^2),3,3) } } cat("Markov chain algorithm ran for ",NUMIT,"\n", "acc.rate for a12d=",acc_a12d/(NUMIT-1),"\n","acc.rate for a3=",acc_a3/(NUMIT-1),"\n", "acc.rate for a4=",acc_a4/(NUMIT-1),"\n","acc.rate for a5=",acc_a5/(NUMIT-1),"\n","acc.rate for a6=",acc_a6/(NUMIT-1),"\n","acc.rate for a7=",acc_a7/(NUMIT-1),"\n", "acc.rate for cx=",acc_cx/(NUMIT-1),"\n","acc.rate for cy=",acc_cy/(NUMIT-1),"\n") return(mchain) } /*Function call for MCMC updates*/ a1_hat1<-MHsampler(25000,s=0.002,s1=0.001,s2=0.27,s3=0.003) /* Estimation of autocorrelation parameter(s) K: size of the sub-regions in the image residual: residual image obtained after removing the trend */ bb=seq(1,1024,K) a_entry=bb[-length(bb)]+1 aa=c(1,a_entry) rhob=matrix(0,length(aa),length(bb)) for (i in 1:length(aa)){ for (j in 1:length(bb)){ block=as.vector(residual[aa[i]:bb[i],aa[j]:bb[j]]) if (any(is.na(block))==FALSE) { rhob[i,j]=sp.correlogram(nbd,block,order=5,method="corr",style="W",zero.policy=T)\$res[1] } else { rhob[i,j]=NA } } } /***********************************************/ /* Spatial analysis of simulated image data */ /***********************************************/ /*define function to remove NA's */ new_sum=function(x) { l=sum(x) is.na(l) } /*define nbd of each pixel: excluding the edge pixels*/ x2=2:99 g1=expand.grid(x=x2,y=x2) nbd=matrix(c(g1\$x+100*(g1\$y-1),ifelse(g1\$y>1,100*(g1\$y-2)+g1\$x,0),ifelse(g1\$x>1,100*(g1\$y-1)+g1\$x-1,0),ifelse(g1\$x+1<=100,g1\$x+1+100*(g1\$y-1),0),ifelse(g1\$y+1<=100,g1\$x+100*g1\$y,0)),9604,5) ### Likelihood function for spatial dependent case ### lpostf<-function(a1,a2,a3,a4,a5,a6,a7,a8,y=data,cxx,cyy,d,nd=nbd){ cx=floor(cxx) cy=floor(cyy) del_cx=cxx-cx del_cy=cyy-cy I<-rep(0,10000) res<-rep(0,10000) S<-rep(0,10000) Kfunc<-rep (0,10000) x1<-seq(0.001,sqrt(2)*0.075,length=100) y1<-seq(0.001,sqrt(2)*0.075,length=100) g=expand.grid(x=x1,y=y1) S<-sqrt((g\$x-x1[cx]+0.000001-del_cx*0.001061273)^2+(g\$y-y1[cy]+0.000001-del_cy*0.001061273)^2) a<-c(a1,a2,a3,a4,a5,a6,a7,a8)/100000 Kfunc<-(a[1]*b(S,1,d)+a[2]*b(S,2,d)+a[3]*b(S,3,d)+a[4]*b(S,4,d)+a[5]*b(S,5,d)+a[6]*b(S,6,d)+a[7]*b(S,7,d)+a[8]*b(S,8,d)) U<-ifelse((Kfunc > 0)==T,Kfunc,0.00001) I<-(U*100000)/S res<-(y-I) Vinv=1/as.vector(I) suml=0 Res1=res[nd[,1]] Res2=res[nd[,2]] Res3=res[nd[,3]] Res4=res[nd[,4]] Res5=res[nd[,5]] Var1=Vinv[nd[,1]] Var2=Vinv[nd[,2]] Var3=Vinv[nd[,3]] Var4=Vinv[nd[,4]] Var5=Vinv[nd[,5]] A1=(Var1+0.25*(Var2+Var3+Var4+Var5))*(Res1^2) A2=(Var2+0.015625*Var1)*(Res2^2) A3=(Var3+0.015625*Var1)*(Res3^2) A4=(Var4+0.015625*Var1)*(Res4^2) A5=(Var5+0.015625*Var1)*(Res5^2) A12=2*(0.5*Var2+0.125*Var1)*(Res2*Res1) A13=2*(0.5*Var3+0.125*Var1)*(Res3*Res1) A14=2*(0.5*Var4+0.125*Var1)*(Res4*Res1) A15=2*(0.5*Var5+0.125*Var1)*(Res5*Res1) A23=2*(0.015625*Var1)*(Res2*Res3) A24=2*(0.015625*Var1)*(Res2*Res4) A25=2*(0.015625*Var1)*(Res2*Res5) A34=2*(0.015625*Var1)*(Res3*Res4) A35=2*(0.015625*Var1)*(Res3*Res5) A45=2*(0.015625*Var1)*(Res4*Res5) detV=log(Var1)+log(Var2)+log(Var3)+log(Var4)+log(Var5) llike=0.125*((A1+A2+A3+A4+A5+A12+A13+A14+A15+A23+A24+A25+A34+A35+A45)-0.5*detV) suml=sum(llike,na.rm=T) prior_part<-sum(log(a))+sum((0.22*(log(a[1:8])-1)^2))+2*d-109*log(d)+(((cxx-50)^2)/8)+(((cyy-50)^2)/8) logpost<-(-(prior_part+suml)) return(logpost) } /********************************************/ /* P(r) from estimated coeffcients */ /********************************************/ /* trapezoidal function: to calculate area under curve */ trap.rule <- function(x,y) { sum(diff(x)*(y[-1]+y[-length(y)]))/2 } /* Plot reference P(r): for simulated data */ /* a<-c( 5.281465e-05, 4.178242e-05, 1.686120e-05 ,5.454271e-06, 3.379092e-06, 3.007528e-06, 2.055125e-06, 8.427587e-07) d<-58.5 #Dmax */ r<-seq(0,d,length=500) # evaluate P(r) on 500 discrete r values between (0,d) Pr_ref=8*pi*r*(a[1]*sin((pi*r*1)/d)+a[2]*sin((pi*r*2)/d)+a[3]*sin((pi*r*3)/d)+a[4]*sin((pi*r*4)/d)+a[5]*sin((pi*r*5)/d)+a[6]*sin((pi*r*6)/d)+a[7]*sin((pi*r*7)/d)+a[8]*sin((pi*r*8)/d)) Pr_ref_const=trap.rule(r,Pr_ref) # normalizing constant plot((Pr_ref/Pr_ref_const)~r, xlab="r", ylab="P(r)",col="gray",type="l",lwd=1,ylim=c(0,0.04)) /* Estimated P(r)'s using esimated parameters from MCMC Numcurves: 10000 # of sampled posterior P(r)'s a1_hat1: MCMC samples of parameters */ Pr_est<-matrix(NA,nrow=500,ncol=Numcurves) R<-matrix(NA,nrow=500,ncol=Numcurves) Pr_std<-matrix(NA,nrow=500,ncol=Numcurves) ##standardized P(r) such that the area under the curve sum to 1 for (k in 1:Numcurves){ R[,k]=seq(0,a1_hat[11,k+15000],length=500) Pr_est[,k]=8*pi*R[,k]*(a1_hat1[1,k+15000]*sin((pi*R[,k]*1)/R[500,k])+a1_hat1[2,k+15000]*sin((pi*R[,k]*2)/R[500,k])+a1_hat1[3,k+15000]*sin((pi*R[,k]*3)/R[500,k])+a1_hat1[4,k+15000]*sin((pi*R[,k]*4)/R[500,k])+a1_hat1[5,k+15000]*sin((pi*R[,k]*5)/R[500,k])+a1_hat1[6,k+15000]*sin((pi*R[,k]*6)/R[500,k])+a1_hat1[7,k+15000]*sin((pi*R[,k]*7)/R[500,k])+a1_hat1[8,k+15000]*sin((pi*R[,k]*8)/R[500,k])) Pr_std[,k]=Pr_est[,k]/trap.rule(R[,k],Pr_est[,k]) } /* Integrated error with respect to reference P(r) for both 2-D and 1-D method x=matrix of estimated parameter values from MCMC updates a=coefficients for reference P(r) am=coefficients from 1-D method */ integrated_error=function(x,a,am){ d=a[length(a)] ll=dim(x)[2] dm=am[length(am)] r=seq(0,d,length=500) Porig=8*pi*r*(a[1]*sin((pi*r*1)/d)+a[2]*sin((pi*r*2)/d)+a[3]*sin((pi*r*3)/d)+a[4]*sin((pi*r*4)/d)+a[5]*sin((pi*r*5)/d)+a[6]*sin((pi*r*6)/d)) const_mcmc=trap.rule(r,Porig) #cat("c=",const_mcmc) rm=seq(0,dm,length=500) Pmoore=8*pi*rm*(am[1]*sin((pi*rm*1)/dm)+am[2]*sin((pi*rm*2)/dm)+am[3]*sin((pi*rm*3)/dm)+am[4]*sin((pi*rm*4)/dm)+am[5]*sin((pi*rm*5)/dm)+am[6]*sin((pi*rm*6)/dm)) const_moore=trap.rule(rm,Pmoore) for (k in 1:ll){ R[,k]=seq(0,x[9,k+n],length=500) PR[,k]=8*pi*R[,k]*(x[1,k+n]*sin((pi*R[,k]*1)/R[500,k])+x[2,k+n]*sin((pi*R[,k]*2)/R[500,k])+x[3,k+n]*sin((pi*R[,k]*3)/R[500,k])+x[4,k+n]*sin((pi*R[,k]*4)/R[500,k])+x[5,k+n]*sin((pi*R[,k]*5)/R[500,k])+x[6,k+n]*sin((pi*R[,k]*6)/R[500,k])) num[k]=trap.rule(R[,k],PR[,k]) abs_diff[k]=sum(abs((PR[,k]/num[k])-(Porig/const_mcmc))) sum_moore=sum(abs((Pmoore/const_moore)-(Porig/const_mcmc))) } return(list(abs_diff=abs_diff, sum_moore=sum_moore)) } /* call integrated_error */ R=matrix(0,500,dim(a1_hat1)[2]) #a1_hat1: matrix of estimated parameter values from 2-D method PR=matrix(0,500,dim(a1_hat1)[2]) PR1=matrix(0,500,dim(a1_hat1)[2]) abs_diff=rep(0,dim(a1_hat1)[2]) num=rep(0,dim(a1_hat1)[2]) amoore: #coefficients from 1D method error_2D=integrated_error(a1_hat1,a,amoore)\$abs_diff error_1D=integrated_error(a1_hat1,a,amoore)\$sum_moore /* Generate crossvalidation dataset (mask) for experimental image */ num=133 # of cells sampled in each x- and y- direction x1=sort(sample(1:1024,num)) y1=sort(sample(1:1024,num)) l1=numeric(1024) l2=numeric(1024) l1[x1]=1 a1=l1 l2[y1]=1 a2=l2 g=expand.grid(x=a1,y=a2) sumg=g\$x+g\$y tt=ifelse(sumg >1,NA,0) ##keep all cells except the sampled cells for test test_mask=matrix(tt,1024,1024) vv=ifelse(sumg >1,0,NA) ##keep only the sampled cells for validation valid_mask=matrix(vv,1024,1024) valid_mask=matrix(vv,1024,1024)