KF code for BM type of black dot

 KF code for BM type of black dot with simulation data and read data

Animation_realData.r

@@ -0,0 +1,45 @@
+
+source('Data_real.r')
+
+postate<-read.csv(file="KappaState2019-05-16.txt",sep="",header=FALSE)
+
+end=max(dim(postate) )
+burnin = 7000 #30000
+mstate<-apply(postate[burnin:end,],2,mean)
+
+
+#max( ax[-which( is.na(ax) ) ] )
+#min( ax[-which( is.na(ax) ) ] )
+
+pl_state = mstate
+pl_state[which(pl_state>1.5)] = 17      # BM
+pl_state[which(pl_state<1.5)] = 19      # OU
+#pl_state[which( is.na(ax) ) ] = NA
+
+pl_state = matrix(c(pl_state),ncol=5)
+
+break_seconds = 1.5
+
+#animation_plot<-function( break_seconds)
+#{
+     plot(0,0,xlim=c(-1,57),ylim=c(-1,45),xlab="x-dimension",ylab="y-dimension",col='white')
+
+       for(qq in 1:46)
+      {
+          lines( ax[qq,1],ay[qq,1],type="p",pch= pl_state[qq,1] ,col="red",cex=1.5 )
+          lines( ax[qq,2],ay[qq,2],type="p",pch= pl_state[qq,2],col="blue",cex=1.5  )
+          lines( ax[qq,3],ay[qq,3],type="p",pch= pl_state[qq,3],col="forestgreen",cex=1.5  )
+          lines( ax[qq,4],ay[qq,4],type="p",pch= pl_state[qq,4],col="black" ,cex=1.5 )
+          lines( ax[qq,5],ay[qq,5],type="p",pch= pl_state[qq,5],col="purple" ,cex=1.5 )
+
+       Sys.sleep(break_seconds)
+          lines( ax[qq,1],ay[qq,1],type="p",pch= pl_state[qq,1],col="white",cex=2  )
+          lines( ax[qq,2],ay[qq,2],type="p",pch= pl_state[qq,2],col="white",cex=2  )
+          lines( ax[qq,3],ay[qq,3],type="p",pch= pl_state[qq,3],col="white",cex=2  )
+          lines( ax[qq,4],ay[qq,4],type="p",pch= pl_state[qq,4],col="white",cex=2  )
+          lines( ax[qq,5],ay[qq,5],type="p",pch= pl_state[qq,5],col="white",cex=2  )
+       }
+#}
+
+
+

Data_real.r

@@ -0,0 +1,54 @@
+
+
+
+library(mvtnorm)
+library(Matrix)
+set.seed(3809)
+
+# Read data
+
+raw <- read.csv("regular.scaled.njaarke.csv")  ## datafile
+raw <- raw[-1,-1]
+
+dates <- as.POSIXct(raw[,1], tz = "GMT")  
+raw <- cbind(dates, raw[,-1])
+
+# time difference between obs
+diff_t <- diff(as.numeric(raw[,1]))/ 3600    # time in mins. Use /3600 to use diff in hours
+dt=diff_t[1]
+
+xy <- as.matrix(raw[,-1])
+nai <- ncol(xy)/2
+Nsamp = 
+alen <- nrow(xy)
+
+animaly<-xy[,2*(1:nai)] / 20          # All y coordinates for each individual
+animalx<-xy[,2*(1:nai)-1] / 20        # All x coordinates for each individual 
+
+# Coordinates for the black dot as an average of all other individuals. 
+black_x <- mean(na.omit(animalx[1,]))
+black_y <- mean(na.omit(animaly[1,])) 
+
+# Initial parameter values
+theta_x <- mean(animalx[1,][!is.na(animalx[1,])])
+theta_y <- mean(animaly[1,][!is.na(animaly[1,])])  
+theta <- c(theta_x,theta_y)
+
+swlamda=c(0.1,0.5)
+
+## rename the data for inferencing code 
+ax = animalx
+ay = animaly
+all_t = seq(0,(alen-1)*2 ,dt)
+all_s = matrix(c(1), nrow = alen, ncol = nai+1 ) ## nai+1 for animal and black dot
+sw_ind = matrix(c("SP"),nrow = alen)
+
+#black = rbind(black_x,black_y)
+#mxsamp = max(all_t)
+
+deltat = dt
+
+
+
+
+

Data_sim.r

@@ -0,0 +1 @@
+library(mvtnorm)library(Matrix)set.seed(3)alpha <- 1.2 #0.8beta  <- 1e-6 rho   <- 5  #5.2sigma <- 0.7   # 2.7Bsigma <- 2swlamda=c(0.1,0.4) #switch rate     attract to unattract       unattract to attract#delata T, the sampling time or the grid time inverval  # Nsamp how many samples we want  deltat=2    dt = deltat ## used for inference code in other files  Nsamp = 50  samp = 1   # how many animalsn = 5#initial black dot and animal locationinit_x = c( 1, 1.5, 1.5, -1, 1.7, 1.7)#, -1.2)init_y = c( 1, 1.5, 1,   -1, 1.7, 1.2)#, -1.2)# initial statesstate<-rep(c(1),n+1)    # state 1 attract to black dot         state 2 unattract to black dot (brownian motion)                        # first element of state is uselessPx = init_xPy = init_yswt_ls <- list( matrix(c(0)),matrix(c(0)),matrix(c(0)),matrix(c(0)),matrix(c(0)) ) # 5 matrix elements for 5 animals  sws_ls <- swt_ls   sw_all <- list(matrix(c(0)), matrix(c(0)),matrix(c(0)))  sampl <- list( matrix(c(0)) )  sw_ind <- c('SP')  all_t <- c(0)  all_s <- state  all_x <- c(Px)  all_y <- c(Py)  swt_old =0                                 # A matrix   coefficient matrix of SDEAid= diag(n+1)A=Aid*(-alpha)A[1:n+1,1] = alphaA[1,1]= -betacovA=rho^2*alpha/(2*beta*(alpha+beta))varA=sigma^2/(2*alpha)+rho^2*alpha/(2*beta*(alpha+beta))lamda =Aid*varA  lamda[1:(n+1),1:(n+1)] =lamda[1:(n+1),1:(n+1)]+covA-covA*Aid[1:(n+1),1:(n+1)]lamda[1,1] = rho^2/(beta*2)## record states at sampling point for all animalstate_sp = array(c(1),dim=c(Nsamp,n))  location_sp = array( (0),dim = c(Nsamp,n+1,2) )#rest T ,how far from current time to the grid time point resT= deltat## propose first switch timeswt=rexp(1,n*swlamda[1])  #switch time  resT = swt - deltatmove2 = swtmove1 = deltatqq = 0tt = 0   repeat{##case 1, switch animal state	if(resT<0)     ## switching happen first , before sampling	{    #cat("switching happen first############# \n")	### move to switch point		t = move2  		#cat("movement at time",t,"\n")		### make the move		ouindex=which(state==1)		bmindex=which(state==2)			    		if(length(ouindex)>1)  		{      ## move black dot and animal with OU			 csigma=lamda[ouindex,ouindex]-expm(A[ouindex,ouindex]*t)%*%lamda[ouindex,ouindex]%*%expm(t(A[ouindex,ouindex])*t)   # var  c- ACA			 csigma=trunc(csigma*10^5)/10^5			 Casigma=as.matrix(csigma)     			     			 mux<-expm(A[ouindex,ouindex]*t)%*%Px[ouindex]                           # mean			 muy<-expm(A[ouindex,ouindex]*t)%*%Py[ouindex]			 cat("animal",ouindex," move OU \n")			   			 Px[ouindex]=rmvnorm(1,as.matrix(mux),Casigma,"chol")			 Py[ouindex]=rmvnorm(1,as.matrix(muy),Casigma,"chol")   		  }		if(length(bmindex)>0)		{     ## move animal with BM			#cat("animal",bmindex,"move BM \n")		   			bmcov=diag(Bsigma^2,length(bmindex))*t			   			Px[bmindex]=rmvnorm(1,Px[bmindex],bmcov,"chol")			Py[bmindex]=rmvnorm(1,Py[bmindex],bmcov,"chol")		}			   			    		## pick the switch animal  		de=sum(swlamda[state])		# SwitchAnimal =sample(c(2:(n+1) ),1,FALSE,prob=c( swlamda[state[2]]/de,swlamda[state[3]]/de,swlamda[state[4]]/de,swlamda[state[5]]/de,swlamda[state[6]]/de,swlamda[state[7]]/de ))    SwitchAnimal =sample( c(2:(n+1) ),1,FALSE,prob=c( swlamda[ state[ 2:(n+1) ] ]/de) )		## switch animal state 		state[SwitchAnimal] = 3 - state[SwitchAnimal]        ## demo state and switch time        swt_ls[[SwitchAnimal-1]] <- cbind( swt_ls[[SwitchAnimal-1]], swt )         sws_ls[[SwitchAnimal-1]] <- cbind( sws_ls[[SwitchAnimal-1]], state[SwitchAnimal] ) 		sw_all[[1]]<- cbind( sw_all[[1]],swt )		sw_all[[2]]<- cbind( sw_all[[2]],state[SwitchAnimal] )		sw_all[[3]]<- cbind( sw_all[[3]], (SwitchAnimal-1) )        cat(tt,"iteration ", "switch", swt ,"\n")                sw_ind= rbind(sw_ind,'SW')        all_s = rbind(all_s,state)        all_t = cbind(all_t,swt+ sum(swt_old) )        all_x = rbind( all_x, rep(NA,n+1) )        all_y = rbind( all_y, rep(NA,n+1) )        swt_old = c(swt_old,swt)		## update index   		ouindex=which(state==1)		bmindex=which(state==2)		## propse next switch 		# cat("ouindex",ouindex,"bmindex",bmindex,"sum","\n")		swt=rexp(1, sum(length(ouindex)*swlamda[1],length(bmindex)*swlamda[2]) ) 		#cat(swt,"next switching time \n")		move1 = abs(resT)		move2 = swt		resT = swt - abs(resT)	#cat(resT,"rest time to next sampling point \n")	#cat("state",state[2]," ",state[3]," ",state[4]," ",state[5]," ",state[6]," ",state[7]," \n \n")	     cat(tt,"iteration ", "switch", swt ,"\n") 	 }###case 2, do not switch animal state	if(resT > 0)            #sampling happan first  resT>0	 {    # cat("sampling happen first############# \n")	### move to sampling point	     t = move1  	# cat("movement at time",t,"\n")	### make move		  ouindex=which(state==1)		  bmindex=which(state==2)		  if(length(ouindex)>1)  		  {		  csigma=lamda[ouindex,ouindex]-expm(A[ouindex,ouindex]*t)%*%lamda[ouindex,ouindex]%*%expm(t(A[ouindex,ouindex])*t)   # var  c- ACA		  csigma=trunc(csigma*10^5)/10^5		  Casigma=as.matrix(csigma)     		     		  mux<-expm(A[ouindex,ouindex]*t)%*%Px[ouindex]                           # mean		  muy<-expm(A[ouindex,ouindex]*t)%*%Py[ouindex]		   #cat("animal",ouindex,"move OU \n")		   		  Px[ouindex]=rmvnorm(1,as.matrix(mux),Casigma,"chol")		  Py[ouindex]=rmvnorm(1,as.matrix(muy),Casigma,"chol")		  }		 if(length(bmindex)>0)		   {		   #cat("animal",bmindex,"move BM \n")		   bmcov=diag(Bsigma^2,length(bmindex))*t		   		   Px[bmindex]=rmvnorm(1,Px[bmindex],bmcov,"chol")		   Py[bmindex]=rmvnorm(1,Py[bmindex],bmcov,"chol")		   }		   	   ### recording location and state information# 	    if(samp == Nsamp)            {              break                } else             {              # save location and state at observation points              location_sp[samp+1,,] = cbind(Px,Py)              state_sp[(samp+1),(1:n)] = state[2:(n+1)]                            # demo state and sampling time              demot = samp*2              cat(tt,"iteration sample at", demot ,"\n")                samp = samp+1              sampl = cbind( sampl, demot)              sw_ind= rbind(sw_ind,'SP')              all_s=  rbind(all_s,state)        	    all_t = cbind(all_t,demot)        	    all_x = rbind( all_x, Px )              all_y = rbind( all_y, Py )               #cat("state",state[2]," ",state[3]," ",state[4]," ",state[5]," ",state[6]," ",state[7]," \n \n")              }	     move2 = abs(resT)	     move1 = deltat	     resT = resT - deltat	 #cat(resT,"rest time to next sampling point \n")	  }  tt=tt+1}ss<- state_sp #t( rbind(s1,s2,s3,s4,s5,s6) )ax <- location_sp[,-1,1] #cbind(a1[1,], a2[1,], a3[1,], a4[1,], a5[1,], a6[1,])ay <- location_sp[,-1,2] #cbind(a1[2,], a2[2,], a3[2,], a4[2,], a5[2,], a6[2,])  animalx = axanimaly = ay theta_x = 0 theta_y = 0 theta=c(theta_x,theta_y)animation_plot<-function( break_seconds){   n = dim(location_sp)[2]-1   colour_l= c('black','red','blue','green','red','blue','green')   pl_state = state_sp   pl_state[which(pl_state>1.5)] = 17      # BM   pl_state[which(pl_state<1.5)] = 19      # OU     plot(0,0,xlim=c( min(location_sp[,,1])-2, max(location_sp[,,2])+2 ),ylim=c( min(location_sp[,,2])-2, max(location_sp[,,2])+2 ),xlab="x-dimension",ylab="y-dimension")       for(qq in 1:Nsamp)      {         for(ia in 1:n)         {          lines(location_sp[qq,ia,1],location_sp[qq,ia,2],type="p",pch=pl_state[qq,ia],col=colour_l[ia])         }       Sys.sleep(break_seconds)                for(ia in 1:n)         {          lines(location_sp[qq,ia,1],location_sp[qq,ia,2],type="p",pch=pl_state[qq,ia],col="white",cex=2)         }      }}