fix lambda problem

mu2013 · Apr 16, 2020 · 7e62374 · 7e62374
1 parent 059830a
commit 7e62374
Show file tree

Hide file tree

Showing 5 changed files with 111 additions and 30 deletions.
diff --git a/GenStates.r b/GenStates.r
@@ -213,7 +213,7 @@ iterkappa = iterkappa +1
 
 ## if we cannot make the proposed state the same as the end time point in 5000 iterations, 
 ## we just reject the proposed states and use the old one. 
-  if(iterkappa == 5000)   
+  if(iterkappa == 3000)   
    {
     iterbreak = 2
     break
@@ -300,6 +300,7 @@ if( iterbreak == 1 )
 }
 
 
+
   # Kappa = max(swlamda)
   # pro_all_time=c(0)
   # repeat
@@ -320,3 +321,14 @@ if( iterbreak == 1 )
   #   }         
   # }
 
+# all_pon = pro_st$pro_pon_ind[ which( pro_st$pro_pon_ind != 'sp') ]
+# all_pon[which(all_pon=='pn')] = 0
+# all_pon[which(all_pon=='OB')] = -1
+# all_pon[which(all_pon=='BO')] = 1
+# all_pon = as.numeric(all_pon)
+# all_pon =cumsum(all_pon)
+
+
+
+
+
diff --git a/Main_Inference.r b/Main_Inference.r
@@ -1,9 +1,10 @@
 
 source('Setup.r')
 source('GenStates.r')
+source('lam_lik.r')
 set.seed(3809)
 
-iter = 50  #  50000
+iter =  50000
 
 nstate = all_s[,-1] 
 nsw_ind = sw_ind
@@ -27,7 +28,17 @@ osp_state = ostate[which(osw_ind=='SP'),]    ## states at sampling points
 opon_ind = pro_st$pro_pon_ind                ## switch indicators at all proposed swiching and sampling points 
 opon_time = pro_st$pro_pon_time              ## time at all proposed swiching and sampling points 
 opon_probAll = pro_st$pro_pon_probAll        ## probability at all proposed swiching and sampling points
-#olnlik_lam = -Inf
+
+ #  sumstate = apply(ostate,1,sum)[-1]
+ #  test = opon_ind[which(opon_ind!='pn')]
+
+ #  test[which(test=='sp')] = 0
+ #  test[which(test=='OB')] = 1
+ #  test[which(test=='BO')] = -1
+ #  test = as.numeric(test)
+
+ # if( sum( sumstate-c(cumsum(test)+5) ) !=0 )
+ # { break }
 
 
 ptm <- proc.time()
@@ -36,20 +47,15 @@ ptm <- proc.time()
 for(j in 1:iter)
 #for(j in 25001:40000)
 {
-## propose new switching rates ##########
-nlambdaOB <-  rnorm(1,swlamda[1],0.05) 
-nlambdaBO <-  rnorm(1,swlamda[2],0.1) 
 
-if( 0.02<nlambdaOB & nlambdaOB <nlambdaBO & nai*nlambdaBO<Kappa)
-  {
-     nswlamda=c(nlambdaOB ,nlambdaBO)
+  for(jj in 1:5){  ## state estimation loop 
  ## propose states list for all animals with switching time
      uplen = 3  #3 is the minimum, in fact we only update the mid point, the start and end states are fixed
      #start_time = sample( seq(2,94,by=2),1)
      start_time = sample( seq(2,(mxsamp - (uplen-1)*deltat),by=deltat ),1)
      end_time = start_time+(uplen-1)*deltat
 
-     pro_st = GenStates( nswlamda,deltat,  uplen,start_time, osp_state,ostate,otime,osw_ind,opon_ind,opon_time)
+     pro_st = GenStates( swlamda,deltat,  uplen,start_time, osp_state,ostate,otime,osw_ind,opon_ind,opon_time)
       nstate = pro_st$pro_state
       ntime = pro_st$pro_time
       nsw_ind = pro_st$pro_sw_ind
@@ -60,60 +66,81 @@ if( 0.02<nlambdaOB & nlambdaOB <nlambdaBO & nai*nlambdaBO<Kappa)
 
     if(iterbreak == 1)
     {
-      nopon_probAll = pro_st$pro_pon_probAll
-      ## assign a constant switching proability, this part of ln likelihood can be cancelled
-      #nlnlik_lam = sum(log( as.numeric(  opon_probAll[ which(opon_ind!='sp') ]   )  ))
-
       old_par = c(alpha,beta,rho,sigma,theta,Bsigma,swlamda)
-      par =  c(alpha,beta,rho,sigma,theta,Bsigma,nswlamda)
+      par =  c(alpha,beta,rho,sigma,theta,Bsigma,swlamda)
 
-      resst = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind, olik, nstate,ntime,nsw_ind, staccept)
+      resst = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind,opon_ind,opon_time,olik, nstate,ntime,nsw_ind, npon_ind,npon_time ,staccept)
       olik = resst$olik
-      #olnlik_lam = resst$olnlik_lam
 
       osp_state = resst$SPstate
-      swlamda = resst$par[8:9]
+      #swlamda = resst$par[8:9]
 
       ostate = resst$ostate
       otime = resst$otime
       osw_ind = resst$osw_ind
+
+      opon_ind = resst$opon_ind
+      opon_time = resst$opon_time
+
       staccept = resst$accept
     } 
-
-     } else{    
-    swlamda = swlamda
   }
 
-#########################################
+## propose new switching rates ##########
+nlambdaOB <-  rnorm(1,swlamda[1],0.05) 
+nlambdaBO <-  rnorm(1,swlamda[2],0.1) 
+nswlamda=c(nlambdaOB ,nlambdaBO)
+
+if( 0.02<nlambdaOB & nlambdaOB <nlambdaBO & nai*nlambdaBO<Kappa)
+  {
+#########################################  switching parameter
+
+      lam_oldlik = lam_lik( swlamda , opon_ind)
+      lam_newlik = lam_lik( nswlamda , opon_ind)
+
+      lam_sigHR <- exp(lam_newlik - lam_oldlik ) 
+        if(runif(1) < lam_sigHR) 
+         {
+           swlamda = nswlamda
+           lamaccept = lamaccept+1
+          } 
+
+#########################################  
+   }
+
+
+#########################################  diffusion parameter
 
   nalpha<-  rnorm(1,alpha,prop.alpha)   # 0.1
   old_par = c(alpha,beta,rho,sigma,theta,Bsigma)
   par = c(nalpha,beta,rho,sigma,theta,Bsigma)
-  resa = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind, olik, ostate,otime,osw_ind, alaccept)
+  resa = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind,opon_ind,opon_time,olik, ostate,otime,osw_ind,npon_ind,npon_time,alaccept)
+
   alpha = resa$par[1]
   olik = resa$olik
   alaccept = resa$accept
 
     nrho<-( rnorm(1,rho,prop.rho) )   
     old_par = c(alpha,beta,rho,sigma,theta,Bsigma)
     par = c(alpha,beta,nrho,sigma,theta,Bsigma)
-    resr = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind, olik, ostate,otime,osw_ind, rhaccept)
+    resr = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind,opon_ind,opon_time,olik, ostate,otime,osw_ind,npon_ind,npon_time,rhaccept)
+
     rho = resr$par[3]
     olik = resr$olik
     rhaccept = resr$accept
 
    nsigma<- ( rnorm(1,sigma,prop.sigma) )   
     old_par = c(alpha,beta,rho,sigma,theta,Bsigma)
     par = c(alpha,beta,rho,nsigma,theta,Bsigma)
-    ress = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind, olik, ostate,otime,osw_ind, sigaccept)
+    ress = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind,opon_ind,opon_time,olik, ostate,otime,osw_ind,npon_ind,npon_time,sigaccept)
     sigma = ress$par[4]
     olik = ress$olik
     sigaccept = ress$accept 
 
    nBsigma<- ( rnorm(1,Bsigma,prop.Bsigma) )  
    old_par = c(alpha,beta,rho,sigma,theta,Bsigma)
    par = c(alpha,beta,rho,sigma,theta,nBsigma)
-   resBs = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind, olik, ostate,otime,osw_ind, Bsigaccept)
+   resBs = Run_KF(par,old_par,osp_state, ostate,otime,osw_ind, opon_ind,opon_time,olik, ostate,otime,osw_ind,npon_ind,npon_time,Bsigaccept)
    Bsigma = resBs$par[7]
    olik = resBs$olik
    Bsigaccept = resBs$accept
@@ -133,7 +160,7 @@ if( 0.02<nlambdaOB & nlambdaOB <nlambdaBO & nai*nlambdaBO<Kappa)
      {
      #cat("iteration",j,"alpha",alaccept,"beta",beaccept,"rho",rhaccept,"sigma",sigaccept,"theta",theaccept,"Bsigma",Bsigaccept,"State", staccept,"\n")
      #cat("iteration",j,"State", staccept,"alpha",alaccept,"beta",beaccept,"\n") # "sigma",sigaccept
-     cat("iteration",j,"State", staccept,"alpha",alaccept,"rho",rhaccept,"Bsigma",Bsigaccept,"sigma",sigaccept,"iterkappa",pro_st$iterkappa,"\n")
+     cat("iteration",j,"State", staccept,"alpha",alaccept,"rho",rhaccept,"Bsigma",Bsigaccept,"sigma",sigaccept,"lambdas",lamaccept,"iterkappa",pro_st$iterkappa,"\n")
      }
 }
 

diff --git a/RunKF.r b/RunKF.r
@@ -1,5 +1,5 @@
 
-Run_KF<-function(par,old_par,osp_state,ostate,otime,osw_ind, olik, nstate,ntime,nsw_ind,acceptm)
+Run_KF<-function(par,old_par,osp_state,ostate,otime,osw_ind,opon_ind, opon_time,olik, nstate,ntime,nsw_ind, npon_ind,npon_time,acceptm)
 {
 #(par,old_par,nstate,osp_state,ostate,ntime,nsw_ind,olik,acceptm)
   nalpha = par[1]
@@ -225,6 +225,8 @@ if(nalpha>threshd&nrho>threshd&nsigma>threshd&nBsigma>threshd)
        ostate = nstate
        otime = ntime
        osw_ind = nsw_ind
+       opon_ind  = npon_ind 
+       opon_time  = npon_time 
    }
    else 
    {
@@ -239,7 +241,7 @@ if(nalpha>threshd&nrho>threshd&nsigma>threshd&nBsigma>threshd)
        sp_state = osp_state
      }
 
-  return( list("par"=c(alpha,beta,rho,sigma,theta,Bsigma,swlamda),"olik"=olik,"accept"=acceptm,"SPstate"=sp_state,"ostate"=ostate,"otime"=otime,"osw_ind"=osw_ind) )
+  return( list("par"=c(alpha,beta,rho,sigma,theta,Bsigma,swlamda),"olik"=olik,"accept"=acceptm,"SPstate"=sp_state,"ostate"=ostate,"otime"=otime,"osw_ind"=osw_ind,"opon_ind"=opon_ind,"opon_time"=opon_time) )
 
 }
 
diff --git a/Setup.r b/Setup.r
@@ -8,7 +8,7 @@ source("SolveSdeBM.r")
 source('RunKF.r')
 
 
-source('Data_sim.r')
+ source('Data_sim.r')
 #source('Data_real.r')
 
 mxsamp = max(all_t)
@@ -63,6 +63,7 @@ osw = -10^4
 
 
 ################## initialise accept number ################################################## 
+lamaccept = 0
 alaccept  = 0
 #beaccept  = 0
 sigaccept = 0

diff --git a/lam_lik.r b/lam_lik.r
@@ -0,0 +1,39 @@
+
+## input
+## lambda_sw switching rate, ou-bm or bm-ou.
+## pro_pon_ind  the indicator vector of states, pn - potential switch, OB - ou to bm, BO - bm to ou.
+
+## output
+## the loglikelihood given the lambdas.
+
+lam_lik<-function( lambda_sw, pro_pon_ind)
+{
+  lambda_BO = lambda_sw[2]
+  lambda_OB = lambda_sw[1]
+
+  p_BO = lambda_BO /Kappa #
+  p_OB = lambda_OB /Kappa # 
+
+  all_pon = pro_pon_ind[ which( pro_pon_ind != 'sp') ]
+  pn_inx = which(all_pon=='pn')
+  OB_inx = which(all_pon=='OB')
+  BO_inx = which(all_pon=='BO')
+
+  num_pon = all_pon
+  num_pon[pn_inx] = 0
+  num_pon[OB_inx] = -1
+  num_pon[BO_inx] = 1
+  num_pon = as.numeric(num_pon)
+
+  cum_pon =cumsum(num_pon) + nai  ## if nai  in cum_pon, all animal doing OU at that time.
+
+  num_pon[OB_inx] = p_OB
+  num_pon[BO_inx] = p_BO
+  num_pon[pn_inx] = 1 - (p_BO*(nai-cum_pon[pn_inx])  + p_OB*(cum_pon[pn_inx]))
+
+  sum(log(num_pon))
+}
+
+
+
+