Unexpected number of observations

[idontgetoutmuch]

I am generating data using

model Lorenz {
  const rho = 45.92
  const beta = 4.0
  const alpha = 16.0


  state X, Y, Z
  obs X_obs


  sub initial {
    X ~ log_normal(log(1.0), 0.00002)
    Y ~ log_normal(log(1.0), 0.00002)
    Z ~ log_normal(log(1.0), 0.00002)
  }


  sub transition(delta = 0.0001) {
    ode {
      dX/dt = alpha * (Y - X)
      dY/dt = X * (rho - Z) - Y
      dZ/dt = X * Y - beta * Z
    }
  }


  sub observation {
    X_obs ~ normal(X, 0.2)
  }
}

Running the following R gives me what I expect: noisy observations.

library('rbi')
library(ggplot2)


model_file_name <- "LorenzGenerate.bi"


Lorenz <- bi_model(model_file_name)


T <- 0.2
nObs <- 100
init_parameters <- list(X = 1, Y = 1, Z = 1)


synthetic_dataset <- bi_generate_dataset(end_time=T, model=Lorenz,
                                         init=init_parameters,
                                         noutputs = nObs)


synthetic_data <- bi_read(synthetic_dataset)


synthetic_df <- as.data.frame(synthetic_data)


ggplot(synthetic_df, aes(X.time)) +
    geom_path(aes(y = X.value, colour="X.value")) +
    geom_path(aes(y = X_obs.value, colour="X_obs.value")) +
    theme(legend.position="bottom") +
    ggtitle("Lorenz") +
    theme(plot.title = element_text(hjust = 0.5)) +
    xlab("Time") +
    ylab("Value")

Notice that I have 101 observations.

> dim(synthetic_df)

[1] 101   9

Now I infer the state using the following

model LorenzState {

  const rho = 45.92

  const beta = 4

  const h         = 0.1;    // time step

  const delta_abs = 1.0e-3; // absolute error tolerance

  const delta_rel = 1.0e-6; // relative error tolerance

  state X, Y, Z

  state ln_alpha

  param mu, sigma

  noise w

  obs X_obs

  sub parameter {

    mu ~ uniform(12.0, 20.0)

    sigma ~ uniform(0.0, 0.5)

  }

  sub proposal_parameter {

     mu ~ truncated_gaussian(mu, 0.02, 12.0, 20.0);

     sigma ~ truncated_gaussian(sigma, 0.01, 0.0, 0.5);

   }

  sub initial {

    X ~ log_normal(log(1.0), 0.2)

    Y ~ log_normal(log(1.0), 0.2)

    Z ~ log_normal(log(1.0), 0.2)

    ln_alpha ~ gaussian(log(mu), sigma)

  }

  sub transition(delta = h) {

    w ~ normal (0.0, sqrt(h))

    ode(h = h, atoler = delta_abs, rtoler = delta_rel, alg =  'RK4(3)') {

      dX/dt = exp(ln_alpha) * (Y - X)

      dY/dt = X * (rho - Z) - Y

      dZ/dt = X * Y - beta * Z

      dln_alpha/dt = -sigma * sigma / 2 - sigma * w / h

    }

  }

  sub observation {

    X_obs ~ log_normal(X, 0.2)

  }

}

I use the following R

model_file_name_state <- "LorenzState.bi"


LorenzState <- bi_model(model_file_name_state)


bi_state <- libbi(model=LorenzState)
bi_state


bi_state <- filter(bi_state, nparticles = 4, nthreads = 1, end_time = T, noutputs = nObs, obs = synthetic_dataset, init = init_parameters)


bi_file_summary(bi_state$output_file_name)
bi_state
summary(bi_state)

But somehow rather than having 101 outputs

> nObs

[1] 100

I have 138

> bi_file_summary(bi_state$output_file_name)

File /private/var/folders/h1/bkwn2ct12hvb6__dzrk5gwx80000gn/T/Rtmp16Twxv/LorenzStatec4966c8f81a/LorenzState_outputc49f4bfe3d.nc (NC_FORMAT_NETCDF4):

     12 variables (excluding dimension variables):

        double time[nr]   (Contiguous storage)  

        double w[np,nr]   (Contiguous storage)  

        double X[np,nr]   (Contiguous storage)  

        double Y[np,nr]   (Contiguous storage)  

        double Z[np,nr]   (Contiguous storage)  

        double ln_alpha[np,nr]   (Contiguous storage)  

        double mu[]   (Contiguous storage)  

        double sigma[]   (Contiguous storage)  

        8 byte int clock[]   (Contiguous storage)  

        int ancestor[np,nr]   (Contiguous storage)  

        double logweight[np,nr]   (Contiguous storage)  

        double loglikelihood[]   (Contiguous storage)  

     2 dimensions:

        nr  Size:138

        np  Size:4

    3 global attributes:

        libbi_schema: ParticleFilter

        libbi_schema_version: 1

        libbi_version: 1.4.1

Why is this?

In theory this shouldn't be a problem (it looks like the solver has been stopped and restarted a very small time step later)

> output$logweight$time[9] - output$logweight$time[10]

[1] -3.469447e-18

But then xtabs doesn't seem to cope with it

> output <- bi_read(bi_state)

> logw <- xtabs(value ~ time + np, data = output$logweight)

> dim(logw)

[1] 101   4

How do I get the expected number of observations? Alternatively how do I get xtabs to behave itself? I should probably ask the latter on an R forum.
 

[Lawrence Murray]

When doing inference, LibBi will output at the time of each observation, and on a regular grid as specified by noutputs, eliminating the intersection between these two sets.

It seems here that the computations for the grid points are slightly different (as in, the tiniest numerical deviation!) from the observation times that appear in the observation file, and so you get double outputs at some times essentially.

As you say it may not be a problem, just a nuisance. I'm not sure I have a quick fix for you. Postprocess the output perhaps? Use a power of two noutputs to avoid decimal rounding?

Either this or you're somehow passing a different noutputs the second time, e.g. noutputs=100 to generate the data set (indeed this will produce 101 obs by the way, one at time zero and 100 equispaced after), and then noutputs=101 the second time.

Lawrence

--
You received this message because you are subscribed to the Google Groups "LibBi Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to libbi-users...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

När du har kontakt med oss på Uppsala universitet med e-post så innebär det att vi behandlar dina personuppgifter. För att läsa mer om hur vi gör det kan du läsa här: http://www.uu.se/om-uu/dataskydd-personuppgifter/

E-mailing Uppsala University means that we will process your personal data. For more information on how this is performed, please read here: http://www.uu.se/om-uu/dataskydd-personuppgifter/

[idontgetoutmuch]

Thanks for this - the easy solution was to remove noutputs for inference and just get results at observations times