logweights become NaNs

[idontgetoutmuch]

I am trying to estimate a parameter for a Lorenz system. Although this system is well known to be chaotic, it is well behaved at the beginning of its evolution. For example here are 32 paths for X, Y and Z starting at (1.0, 1.0, 1.0) perturbed by N(0, 0.2).

I generate some noisy 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)

  }

}

And then try and infer the parameter alpha by making it a state variable

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)

  }

}

But I am getting NaNs for the logweights

Here's my R

library('rbi')

library(ggplot2)

model_file_name <- "LorenzGenerate.bi"

Lorenz <- bi_model(model_file_name)

T <- 2.0

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")

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, obs = synthetic_dataset, init = init_parameters, ess_rel=1)

And here are the logweights getting smaller and smaller until we get NaNs.

> output$logweight$value[1:20]

 [1]     -8.195611    -28.326456    -44.079191    -56.012739    -84.718022

 [6]   -138.593273   -263.272262   -551.239571  -1135.085810  -2370.284954

[11]  -4742.979880  -8895.140169 -15225.553036 -23177.856321 -31058.553846

[16] -37130.818336 -41081.672274           NaN           NaN           NaN

> output$logweight$value[102:121]

 [1]     -8.195611    -29.026216    -46.347489    -55.193830    -81.421323

 [6]   -138.593273   -261.459513   -537.840374  -1135.085810  -2370.284954

[11]  -4742.979880  -8962.658035 -15225.553036 -23177.856321 -31058.553846

[16] -37130.818336 -41082.880530           NaN           NaN           NaN

> output$logweight$value[203:222]

 [1]     -8.195611    -29.915752    -43.809861    -55.193830    -81.421323

 [6]   -138.593273   -263.271111   -537.840374  -1135.085810  -2370.284954

[11]  -4742.979880  -8896.327921 -15230.638953 -23177.856321 -31058.553846

[16] -37130.818336 -41067.499170           NaN           NaN           NaN

> output$logweight$value[304:323]

 [1]     -8.195611    -28.292957    -43.809861    -55.193830    -81.421323

 [6]   -138.593273   -264.276304   -537.840374  -1135.085810  -2370.284954

[11]  -4742.979880  -8928.661460 -15225.553036 -23177.856321 -31058.553846

[16] -37130.818336 -41006.349718           NaN           NaN           NaN

So it seems that all particles are a long way from the observations and get further away as time goes on. I think I will try running SIR on this model "by hand" so I can see what is going on but if anyone has any advice, please let me know.

[Sebastian Funk]

You've got a different observation model in the two files (normal vs. lognormal). If you use the 'filter' command with 'sample_obs=TRUE' and look at the sampled observations, you'll find that they diverge quickly where X gets too large. If you replace the 'log_normal' with a 'normal' in the second model, you should get better results.

idontgetoutmuch wrote:
> I am trying to estimate a parameter for a Lorenz system. Although this
> system is well known to be chaotic, it is well behaved at the beginning of
> its evolution. For example here are 32 paths for X, Y and Z starting at
> (1.0, 1.0, 1.0) perturbed by N(0, 0.2).
>
>
>

> <https://lh3.googleusercontent.com/-xN23P3cQrfk/WzHpX-7F3jI/AAAAAAAABQ4/9B4WcIE4PJYpgvtSSXgTSHOLoCjsW4JgQCLcBGAs/s1600/LogChart.png>

>
>
> I generate some noisy 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)
> }
> }
>
>

> <https://lh3.googleusercontent.com/-DJjfBbvZ45w/WzHwKwis5UI/AAAAAAAABRE/Dfu_pmC-T5IS8MYxA_YmF-gxrQp4nSmvQCLcBGAs/s1600/LorenzNoisy.png>

[idontgetoutmuch]

Great catch 😀 Thanks that has saved me a lot of work.

Now you can see that the X state value is being tracked quite nicely even beyond where the system becomes chaotic

and more importantly the parameter (which I modelled as Brownian motion) is giving a nice and similar answer for both systems even thought the states have diverged.

[idontgetoutmuch]

For reasons I don't understand this chart got missed from my reply