Rprintf in Csnippet - Githubissues

I'm trying to estimate parameters in a multi-vairate epidemic model via the mif2 algorithm. Unfortunately there is often an error message that the dmeasure returns non-finite values during the mif2 process. In your documentation it says that using Rprintf in the dmeasure Csnippet can help debug the error. If I use the Rprintf statement in my Csnippet code there is always an error message during the creation of the pomp-object. Is there an example how to use the Rprintf statement in the Csnippet code?

What is the error message you are getting when you try to build the pomp object with the Rprintf statement in it? It would be most helpful if you could provide a simple example that produces the error.

oh, sorry about that. the dmeasure Csnippet looks like:

dmeasure <- Csnippet("
   lik=dnorm(S_share,N2,rho_S,1)+dnorm(I_share,N3,rho_I,1)+
           dnorm(P_share,N4,rho_P,1)+dnorm(D_share,N5, rho_D,1);
   lik = (log) ? lik : exp(lik);
   Rprintf(lik);
")

...so obviously I'd like to get the value of the likelihood and the error message I get is:

error: passing 'double' to parameter of incompatible type 'const char *' Rprintf(lik);

Very good, thanks. Try

dmeasure <- Csnippet("
             lik=dnorm(S_share,N2,rho_S,1)+
                    dnorm(I_share,N3,rho_I,1)+
                    dnorm(P_share,N4,rho_P,1)+
                    dnorm(D_share,N5, rho_D,1);
             lik = (log) ? lik : exp(lik);
             if (!R_FINITE(lik)) 
                 Rprintf(\"%lg %lg %lg %lg %lg\\n\",rho_S,rho_I,rho_P,rho_D,lik);
")

Note the escape characters in the Rprintf format argument are needed because the whole snippet is a text string.

The thought here is that, by printing the values of the parameters associated with non-finite values of the likelihood, you'll get some insight into what is going wrong. My guess is that your standard deviations are going negative. If this turns out to be right, you should consider using parameter transformation to constrain them to be positive.

For more information on Rprintf and other components of R's C interface, check out the corresponding section in the R extensions manual.

Thanks a lot. It worked. The issue is that some states go beyond any limit. But thanks to you I know where the issue is, thanks.

Glad to have been able to help!

Hey, I wasn't sure if I need to open a new topic or continue this one. Thanks to you I was able to get all the states and likelihoods in every time step due to the Rprintf command. At some point all my particles go zero which means some of my states go beyond any limit after that and an error message occurs that dmeasure returns non-finite values. So I read I can choose a tolerance limit via "tol", but it seems the mif2 algorithm uses the particles with the likelihood below tol anyway. The mif2 command looks like:

 m1 <- mif2(
      test,
      Nmif=50,
      start=theta.guess,
      tol = 1e-10,
      transform=TRUE,
      var.factor=20,
      rw.sd=rw,
      cooling.fraction=0.95,
      Np=10000
    )

at some point the likelikood of one of the particles is: 1.88697e-71 and in the next time step: 2.76405e-143

Not sure what you mean by "all particles go to zero". What becomes zero? If the likelihood of the particles is becoming very small, that is because those particles are incompatible with the data. This, in itself, is not a cause for concern. What happens when you run 'pfilter'?

I mean the likelihood of all particles becomes zero at some point and after that the states start to grow to (+/-) infinity so that the algorithm stops with an error. That means I'm not able to run pfilter with mle-parameter values. With my initial parameter guess the pfilter works fine

I am confused, because if what you claim is correct, then the likelihood at your initial parameter guess is higher than the parameters you get after running mif2. If this is really happening (I note that you don't say whether you've actually tried running pfilter at your initial guess), then it is almost certainly a bug either in your code or in mine.

If in pfilter all particles have likelihoods below tol (which is by default 1e-17), then a warning about "filtering failure" is issued. Are you seeing this warning when running mif2 and pfilter? If not, then you are wrong about all your particles having very low likelihood.

At this point, it would be most helpful were you to provide a self-contained example. Failing that, a complete listing of your code and output would be good.

Yes, I'm seeing the warning about "filtering failure". The thing is if I'm running mif2 there is the error message "Error : ‘mif2.pfilter’ error: ‘dmeasure’ returns non-finite value" and if I'm running mif2 a second time with the same initial guess mif2 seems sometimes to work with just the "filtering failure" warning message mentioned above!?

If I'm running pfilter the likelihood with my initial guess is higher than the likelihood with the parameters after running mif2!?

So as you said I provide you with my code:

Prozess<-Csnippet("
                  // Define variables
                  double M[6];
                  double lambda, nu_S, alpha, nu_I, gamma;
                  double epsilon1, dW, epsilon2, epsilon3, epsilon4, epsilon5;

                  // random values
                  epsilon1=rnorm(0,tau2);
                  epsilon2=rnorm(0,tau3);
                  epsilon3=rnorm(0,tau4);
                  epsilon4=rnorm(0,tau5);
                  epsilon5=rnorm(0,tau6);
                  dW=rgamma((dt/(sigma*sigma)),1/(sigma*sigma));

                  // transition rates
                  lambda=(N3+N5)*exp(a1*unempl_real+a2*HousePriceCounty+a3*spread+epsilon2);                              // Bakterienkonzentration in der Bevölkerung
                  nu_S=b0*exp(b1*unempl_real+b2*HousePriceCountyDelta+b3*spread+epsilon1);                                         // Übertragung S->P
                  nu_I=c0*exp(c1*unempl_real+c2*HousePriceCountyDelta+c3*spread+epsilon3);                                         // Übertragung I->P
                  alpha=d0*exp(d1*unempl_real+d2*HousePriceCountyDelta+d3*spread+epsilon4);                                        // Übertragung I->D
                  gamma=e0*exp(e1*unempl_real+e2*HousePriceCountyDelta+e3*spread+e4*X_1y_ARM+e5*X_30y_FRM+e6*FICO+epsilon5);          // Übertragung I->S 

                  // process model
                  M[0]=N0+((lambda-N0)*(2/tau)*dW)*dt;
                  M[1]=N1+((M[0]-N1)*2/tau)*dt;
                  M[2]=N2+(-bet*N1*N2+gamma*N3-nu_S*N2)*dt;
                  M[3]=N3+(bet*N1*N2-gamma*N3-alpha*N3-nu_I*N3)*dt;
                  M[4]=N4+(nu_S*N2+nu_I*N3)*dt;
                  M[5]=N5+(alpha*N3)*dt;
                  N0=M[0];
                  N1=M[1];
                  N2=M[2];
                  N3=M[3];
                  N4=M[4];
                  N5=M[5];

                  ")
rmeasure <- Csnippet("
                     S_share=rnorm(N2, rho_S);
                     I_share=rnorm(N3, rho_I);
                     P_share=rnorm(N4, rho_P);
                     D_share=rnorm(N5, rho_D);
                     ")
dmeasure <- Csnippet("

                     lik = dnorm(S_share,N2, rho_S,1)+dnorm(I_share,N3, rho_I,1)+dnorm(P_share,N4, rho_P,1)+dnorm(D_share,N5, rho_D,1);
                     lik = give_log ? lik : exp(lik);
                      ")
logtrans<-Csnippet("
                   Trho_S=log(rho_S);
                   Trho_I=log(rho_I);
                   Trho_P=log(rho_P);
                   Trho_D=log(rho_D);
                   Ttau2=log(tau2);
                   Ttau3=log(tau3);
                   Ttau4=log(tau4);
                   Ttau5=log(tau5);
                   Ttau6=log(tau6);
                   Tb0=log(b0);
                   Tc0=log(c0);
                   Td0=log(d0);
                   Te0=log(e0);
                   ")
exptrans<-Csnippet("  
                   Trho_S=exp(rho_S);
                   Trho_I=exp(rho_I);
                   Trho_P=exp(rho_P);
                   Trho_D=exp(rho_D);
                   Ttau2=exp(tau2);
                   Ttau3=exp(tau3);
                   Ttau4=exp(tau4);
                   Ttau5=exp(tau5);
                   Ttau6=exp(tau6);
                   Tb0=exp(b0);
                   Tc0=exp(c0);
                   Td0=exp(d0);
                   Te0=exp(e0);
                   ")

# -----starting values of the states------------------------------
init=c(N0.0=0,N1.0=0,N2.0=FIPS_einzeln[1,"S_share"],N3.0=FIPS_einzeln[1,"I_share"],
N5.0=FIPS_einzeln[1,"D_share"],N4.0=FIPS_einzeln[1,"P_share"])

# -----name of the parameters-------------------------------------------
estpars<-c("rho_S","rho_I","rho_P","rho_D","bet","a1","a2","a3",
"b0","b1","b2","b3","c0","c1","c2","c3",
"d0","d1","d2","d3","e0","e1","e2","e3","e4","e5","e6",
"sigma","tau","tau2","tau3","tau4","tau5","tau6") 

 estpars_deftime<-c("a1","a2","a3","b0","b1","b2","b3","c0","c1","c2","c3",
"d0","d1","d2","d3","e0","e1","e2","e3","e4","e5","e6")  

estpars_std<-c("sigma","tau2","tau3","tau4","tau5","tau6", "tau") 

# -----define "pomp"-object-------------------------------------
test <- pomp(data=FIPS_einzeln,                                # Datentabelle
             times="time",                                     # Spalte mit Zeiteinheiten
             t0=with(FIPS_einzeln,time[1]),                    # Startpunkt des betrachteten Zeitfensters
             rmeasure=rmeasure,                                # Messmodell (versch. Var für jeden State) -> alternativ "alle": rmeasure_all oder "ohne B": rmeasure_oB
             dmeasure=dmeasure,                                # Messmodell (versch. Var für jeden State) -> alternativ "alle": dmeasure_all oder "ohne B": dmeasure_oB
             rprocess=discrete.time.sim(Prozess,delta.t=1),    # Simulation der Dynamiken -> alternativ euler.sim(Prozess,delta.t=1/365.25)
             #rprocess=euler.sim(Prozess,delta.t=1/365.25),             
             statenames=c("N0","N1","N2","N3","N4","N5"),      # Namen der States (N0 nur für "B"-Berechnung nötig?!)
             paramnames=estpars,                               # Parameternamen
             covar=covar,                                      # Tabelle mit Kovariaten
             tcovar="time",                                    # Spalte mit Zeiteinheiten für Kovariaten
             toEstimationScale=logtrans,                       # Transformation der Parameter (bringen wie oben definiert keinen Mehrwert)
             fromEstimationScale=exptrans,                     # Transformation der Parameter (bringen wie oben definiert keinen Mehrwert)
             #initializer=init                                 # Startwerte der States
             params=init
)

# ----- random walk for iterated filtering algorithm -------------
rw<-rep(0.1,length(estpars))
names(rw)<-estpars

# ----- initial parameter guess -------------
    theta.guess<-init
    theta.guess[c("N0.0","N1.0")]<-runif(n=2,min=0,max=0.1)
    theta.guess[c("bet")]<-runif(n=1,min=-0.1,max=0.1)
    theta.guess[estpars_deftime]<-runif(n=length(estpars_deftime),min=-1,max=1)
    theta.guess[estpars_std]<-runif(n=length(estpars_std),min=0,max=1)
    theta.guess["sigma"]<-runif(n=1,min=2,max=4)
    theta.guess["tau"]<-runif(n=1,min=0,max=10)
    theta.guess[c("b0","c0","d0","e0")]<-runif(n=4,min=0,max=1)
    theta.guess[c("rho_S","rho_I","rho_P","rho_D")]<-runif(n=4,min=0.5,max=1)

m1 <- mif2(
      test,
      Nmif=50,
      start=theta.guess,
      tol = 1e-17,
       transform=TRUE,
      var.factor=20,
      rw.sd=rw,
      cooling.fraction=0.95,
      Np=10000
    )

output from first try :

Fehler: in ‘mif2’: particle-filter error:Error : ‘mif2.pfilter’ error: ‘dmeasure’ returns non-finite value

output from second try:

Warnmeldungen:
1: 1 filtering failure occurred in ‘mif2.pfilter’ 
2: 3 filtering failures occurred in ‘mif2.pfilter’ 
3: 3 filtering failures occurred in ‘mif2.pfilter’ 
4: 1 filtering failure occurred in ‘mif2.pfilter’ 
5: 3 filtering failures occurred in ‘mif2.pfilter’ 
6: 3 filtering failures occurred in ‘mif2.pfilter’ 
7: 3 filtering failures occurred in ‘mif2.pfilter’ 
8: 4 filtering failures occurred in ‘mif2.pfilter’ 
9: 12 filtering failures occurred in ‘mif2.pfilter’ 
10: 2 filtering failures occurred in ‘mif2.pfilter’

p_filter<-replicate(n=50,logLik(pfilter(test,param=theta.guess,Np=10000,tol = 1e-17)))

-> output around -614

p_filter2<-replicate(n=50,logLik(pfilter(m1,Np=10000,tol = 1e-17)))

-> output between -4400 and -6800

Thanks a lot for your help!

I wasn't able to run your code. I get error

object 'FIPS_einzeln' not found

Because I can't run it, I can only guess as to what the problems might be. Just looking over your code, I note a number of elements that may not be what you intend.

You attempt to estimate 34 parameters (length(estpars)), yet you only transform 13 of these. If, in its search of parameter space, mif2 proposes negative values for some of these parameters, it may violate implicit assumptions of your model. In particular, it may lead to non-finite likelihoods.
I note that you have what looks like a white-noise term dW=rgamma((dt/(sigma*sigma)),1/(sigma*sigma)) in your process model. This will have expected value dt/sigma^4 and variance dt/sigma^6, which I suspect is not your intention. In particular, a Gamma white noise process with mean dt and variance sigma^2*dt can be had by using dW = rgamma(dt/sigma/sigma,sigma*sigma). For convenience, the package provides the equivalent rgammawn command: dW=rgammawn(sigma,dt).
Again, in the process model Csnippet, you have the line M[0]=N0+((lambda-N0)*(2/tau)*dW)*dt. The increment here has expectation proportional to dt^2, which again is not what you intend, I suspect.

Whether the foregoing explain the problem, I can't know without being able to run the codes.

Ok, so I provide you with the data matrix 'FIPS_einzeln' and the covariate matrix 'covar' in the attachment. I changed the transformation so that the only parameters allowed negative values are these inside the exp-function and I adjusted the process model according to your comments:

Prozess<-Csnippet("
// Define variables
                  double M[6];
                  double lambda, nu_S, alpha, nu_I, gamma;
                  double epsilon1, dW, epsilon2, epsilon3, epsilon4, epsilon5;

                  // random values
                  epsilon1=rnorm(0,tau2);
                  epsilon2=rnorm(0,tau3);
                  epsilon3=rnorm(0,tau4);
                  epsilon4=rnorm(0,tau5);
                  epsilon5=rnorm(0,tau6);
                  dW=rgamma((dt/sigma/sigma),sigma*sigma);

                  // transition rates
                  lambda=(N3+N5)*exp(a1*unempl_real+a2*HousePriceCounty+a3*spread+epsilon2);                              // Bakterienkonzentration in der Bevölkerung
                  nu_S=b0*exp(b1*unempl_real+b2*HousePriceCountyDelta+b3*spread+epsilon1);                                         // Übertragung S->P
                  nu_I=c0*exp(c1*unempl_real+c2*HousePriceCountyDelta+c3*spread+epsilon3);                                         // Übertragung I->P
                  alpha=d0*exp(d1*unempl_real+d2*HousePriceCountyDelta+d3*spread+epsilon4);                                        // Übertragung I->D
                  gamma=e0*exp(e1*unempl_real+e2*HousePriceCountyDelta+e3*spread+e4*X_1y_ARM+e5*X_30y_FRM+e6*FICO+epsilon5);          // Übertragung I->S 

                  // process model
                  M[0]=N0+((lambda-N0)*(2/tau)*dW)*dt;
                  M[1]=N1+((M[0]-N1)*2/tau)*dt;
                  M[2]=N2+(-bet*N1*N2+gamma*N3-nu_S*N2)*dt;
                  M[3]=N3+(bet*N1*N2-gamma*N3-alpha*N3-nu_I*N3)*dt;
                  M[4]=N4+(nu_S*N2+nu_I*N3)*dt;
                  M[5]=N5+(alpha*N3)*dt;
                  N0=M[0];
                  N1=M[1];
                  N2=M[2];
                  N3=M[3];
                  N4=M[4];
                  N5=M[5];     
                  ")
logtrans<-Csnippet("
                   Tbet=log(bet);
                   Trho_S=log(rho_S);
                   Trho_I=log(rho_I);
                   Trho_P=log(rho_P);
                   Trho_D=log(rho_D);
                   Ttau=log(tau);
                   Tsigma=log(sigma);
                   Ttau2=log(tau2);
                   Ttau3=log(tau3);
                   Ttau4=log(tau4);
                   Ttau5=log(tau5);
                   Ttau6=log(tau6);
                   Tb0=log(b0);
                   Tc0=log(c0);
                   Td0=log(d0);
                   Te0=log(e0);
                   ")
exptrans<-Csnippet("
                   Tbet=exp(bet);
                   Trho_S=exp(rho_S);
                   Trho_I=exp(rho_I);
                   Trho_P=exp(rho_P);
                   Trho_D=exp(rho_D);
                   Ttau=exp(tau);
                   Tsigma=exp(sigma);
                   Ttau2=exp(tau2);
                   Ttau3=exp(tau3);
                   Ttau4=exp(tau4);
                   Ttau5=exp(tau5);
                   Ttau6=exp(tau6);
                   Tb0=exp(b0);
                   Tc0=exp(c0);
                   Td0=exp(d0);
                   Te0=exp(e0);
                   ")

The issues still exist. If I'm running mif2 the 'dmeasure' returns non-finite value but the pfilter works with the same initial parameter guess!? The mif2 sometimes stops with the error message: ‘mif2.pfilter’ error: ‘dmeasure’ returns non-finite value and if I'm running it again, the algorithm works fine!?

Covar.xlsx FIPS_einzeln.xlsx

I believe I see what the difficulty was. There were several related issues.

In the following, I have modified your process model line M[0]=N0+((lambda-N0)*(2/tau)*dW)*dt to M[0]=N0+(lambda-N0)*(2/tau)*dW. The former is proportional to dt^2; the latter, to dt.

Two aspects of your codes are in contradiction. You log transform the bet parameter, yet in your initial guess, you sometimes assign it negative values. I am not sure whether your model includes the constraint bet>0 or not. In the following, I have retained this constraint and eliminated the possibility that bet is initially assigned a negative value.

Having addressed these issues, I nevertheless was able to see the phenomenon you were describing. The estimated likelihood at the initial guess was often actually higher than that at the ostensive MLE returned by mif2. To see what was going on, I plotted the MIF2 likelihood (as returned by conv.rec or plot on m1) and observed that it was increasing. I interpreted the large discrepancy between the MIF2 likelihood and the actual likelihood (as estimated by pfilter) as evidence that the likelihood surface for the model with perturbed parameters is quite different from that without. (See FAQ 5.1 for more on this.) Accordingly, I reduced the magnitude of the MIF2 random-walk perturbations (via the rw.sd argument of mif2) and was able to increase the likelihood.

This raises the question of what the appropriate size of the random-walk perturbations should be. It should be fairly obvious that the answer must depend strongly on the nature of the model and the data, so that no precise numerical answer is possible. In general, the perturbations should be on roughly the same scale as important features of the local likelihood landscape. If they are much larger (as I believe was the case in the codes you sent me), the important features are obscured by the noise. If they are much smaller, the efficiency of the maximization algorithm will be low. It's probably safe to say, too, that, other things being equal, the larger the parameter space being searched, the smaller the perturbations will have to be.

require(pomp)
require(magrittr)
require(readxl)

read_excel("FIPS_einzeln.xlsx") %>% as.data.frame() -> FIPS_einzeln
read_excel("Covar.xlsx") %>% as.data.frame() -> covar

Prozess<-Csnippet("
// Define variables
                  double M[6];
                  double lambda, nu_S, alpha, nu_I, gamma;
                  double epsilon1, dW, epsilon2, epsilon3, epsilon4, epsilon5;

                  // random values
                  epsilon1=rnorm(0,tau2);
                  epsilon2=rnorm(0,tau3);
                  epsilon3=rnorm(0,tau4);
                  epsilon4=rnorm(0,tau5);
                  epsilon5=rnorm(0,tau6);
                  dW=rgammawn(sigma,dt);

                  // transition rates
                  lambda=(N3+N5)*exp(a1*unempl_real+a2*HousePriceCounty+a3*spread+epsilon2);                              // Bakterienkonzentration in der Bevölkerung
                  nu_S=b0*exp(b1*unempl_real+b2*HousePriceCountyDelta+b3*spread+epsilon1);                                         // Übertragung S->P
                  nu_I=c0*exp(c1*unempl_real+c2*HousePriceCountyDelta+c3*spread+epsilon3);                                         // Übertragung I->P
                  alpha=d0*exp(d1*unempl_real+d2*HousePriceCountyDelta+d3*spread+epsilon4);                                        // Übertragung I->D
                  gamma=e0*exp(e1*unempl_real+e2*HousePriceCountyDelta+e3*spread+e4*X_1y_ARM+e5*X_30y_FRM+e6*FICO+epsilon5);          // Übertragung I->S 

                  // process model
// AAK: the following line has been edited.  NB: E[dW] = dt
                  M[0]=N0+(lambda-N0)*(2/tau)*dW;
                  M[1]=N1+((M[0]-N1)*2/tau)*dt;
                  M[2]=N2+(-bet*N1*N2+gamma*N3-nu_S*N2)*dt;
                  M[3]=N3+(bet*N1*N2-gamma*N3-alpha*N3-nu_I*N3)*dt;
                  M[4]=N4+(nu_S*N2+nu_I*N3)*dt;
                  M[5]=N5+(alpha*N3)*dt;
                  N0=M[0];
                  N1=M[1];
                  N2=M[2];
                  N3=M[3];
                  N4=M[4];
                  N5=M[5];
// AAK: print some information when errors occur
if (!(R_FINITE(N2)&&R_FINITE(N3))) Rprintf(\"%lg %lg %lg %lg %lg %lg %lg %lg\\n\",N1,N2,N3,bet,gamma,alpha,nu_S,nu_I);
                  ")

logtrans<-Csnippet("
                   Tbet=log(bet);
                   Trho_S=log(rho_S);
                   Trho_I=log(rho_I);
                   Trho_P=log(rho_P);
                   Trho_D=log(rho_D);
                   Ttau=log(tau);
                   Tsigma=log(sigma);
                   Ttau2=log(tau2);
                   Ttau3=log(tau3);
                   Ttau4=log(tau4);
                   Ttau5=log(tau5);
                   Ttau6=log(tau6);
                   Tb0=log(b0);
                   Tc0=log(c0);
                   Td0=log(d0);
                   Te0=log(e0);
                   ")

exptrans<-Csnippet("
                   Tbet=exp(bet);
                   Trho_S=exp(rho_S);
                   Trho_I=exp(rho_I);
                   Trho_P=exp(rho_P);
                   Trho_D=exp(rho_D);
                   Ttau=exp(tau);
                   Tsigma=exp(sigma);
                   Ttau2=exp(tau2);
                   Ttau3=exp(tau3);
                   Ttau4=exp(tau4);
                   Ttau5=exp(tau5);
                   Ttau6=exp(tau6);
                   Tb0=exp(b0);
                   Tc0=exp(c0);
                   Td0=exp(d0);
                   Te0=exp(e0);
                   ")

rmeasure <- Csnippet("
                     S_share=rnorm(N2, rho_S);
                     I_share=rnorm(N3, rho_I);
                     P_share=rnorm(N4, rho_P);
                     D_share=rnorm(N5, rho_D);
                     ")

## AAK: 'dmeasure' is modified to print state variables when errors occur
dmeasure <- Csnippet("
                     lik = dnorm(S_share,N2, rho_S,1)+dnorm(I_share,N3, rho_I,1)+dnorm(P_share,N4, rho_P,1)+dnorm(D_share,N5, rho_D,1);
if (!R_FINITE(lik)) Rprintf(\"%lg %lg %lg %lg\\n\",N2,N3,N4,N5);
                     lik = give_log ? lik : exp(lik);
                     ")

# -----starting values of the states------------------------------
init=c(N0.0=0,N1.0=0,N2.0=FIPS_einzeln[1,"S_share"],N3.0=FIPS_einzeln[1,"I_share"],N5.0=FIPS_einzeln[1,"D_share"],N4.0=FIPS_einzeln[1,"P_share"])

# -----name of the parameters-------------------------------------------
estpars<-c("rho_S","rho_I","rho_P","rho_D","bet","a1","a2","a3","b0","b1","b2","b3","c0","c1","c2","c3","d0","d1","d2","d3","e0","e1","e2","e3","e4","e5","e6","sigma","tau","tau2","tau3","tau4","tau5","tau6")
estpars_deftime<-c("a1","a2","a3","b0","b1","b2","b3","c0","c1","c2","c3","d0","d1","d2","d3","e0","e1","e2","e3","e4","e5","e6")
estpars_std<-c("sigma","tau2","tau3","tau4","tau5","tau6", "tau")

# -----define "pomp"-object-------------------------------------
test <- pomp(data=FIPS_einzeln,
             times="time",
             t0=with(FIPS_einzeln,time[1]),
             rmeasure=rmeasure,
             dmeasure=dmeasure,
             rprocess=discrete.time.sim(Prozess,delta.t=1),
             statenames=c("N0","N1","N2","N3","N4","N5"),
             paramnames=estpars,
             covar=covar,
             tcovar="time",
             toEstimationScale=logtrans,
             fromEstimationScale=exptrans,
             params=init
)

try(plot(simulate(test)))  ## AAK: fails because parameters are not fully specified

# ----- initial parameter guess -------------
theta.guess <- init
theta.guess[c("N0.0","N1.0")] <- runif(n=2,min=0,max=0.1)
## AAK: In your model, "bet" can be negative, but it is log-transformed.
## AAK: The following line eliminates the potential for bet < 0.
theta.guess[c("bet")] <- runif(n=1,min=0,max=0.1)
theta.guess[estpars_deftime] <- runif(n=length(estpars_deftime),min=-1,max=1)
theta.guess[estpars_std] <- runif(n=length(estpars_std),min=0,max=1)
theta.guess["sigma"] <- runif(n=1,min=2,max=4)
theta.guess["tau"] <- runif(n=1,min=0,max=10)
theta.guess[c("b0","c0","d0","e0")] <- runif(n=4,min=0,max=1)
theta.guess[c("rho_S","rho_I","rho_P","rho_D")] <- runif(n=4,min=0.5,max=1)

replicate(5,logLik(pfilter(test,Np=1000,params=theta.guess))) %>% logmeanexp(se=TRUE)
## AAK: returns c. -500

# ----- random walk for iterated filtering algorithm -------------
## AAK: note this is 1/100 its former value
rw <- rep(0.001,length(estpars))
names(rw) <- estpars

m1 <- mif2(
  test,
  Nmif=50,
  start=theta.guess,
  transform=TRUE,
  rw.sd=rw,
  cooling.fraction=0.5,
  cooling.type='geometric',
  Np=1000
)

replicate(5,logLik(pfilter(m1,Np=1000))) %>% logmeanexp(se=TRUE)
## AAK: returns c. -98

require(ggplot2)
require(reshape2)

m1 %>%
    conv.rec() %>%
    melt() %>%
    ggplot(aes(x=Var1,y=value))+
    geom_line()+facet_wrap(~Var2,scales="free_y")
## AAK: see attached plot

issue13b

I hope this helps. Note that while the above certainly increases the likelihood above that of the initial guess, there is more work to be done to find the MLE.

BTW, I hope you will take care in future to format your markdown so as to make it easier for me (and others) to extract the codes.

Thanks a lot for your help!

kingaa / pomp

Rprintf in Csnippet #13