alpha = 16.0rho = 45.92beta = 4
dzdt :: Double -> [Double] -> [Double]dzdt _t [x, y, z] = [ alpha * (y - x) , x * (rho - z) - y , x * y - beta * z ]
m = odeSolve dzdt [1.0, 1.0, 1.0] (LA.vector $ take 1000 [1.0e-4, 2.0e-4..])
*Main> m(1000><3) [ 1.0, 1.0, 1.0 , 1.0000035116198551, 1.0043918006561992, 0.9997002796734803 , 1.0000140385819727, 1.0087832236309844, 0.9994011189915194 , 1.000031569092141, 1.0131743003919085, 0.9991025184087792 , 1.0000560914252676, 1.0175650623511352, 0.99880447839092 , 1.0000875939251799, 1.0219555408657492, 0.9985069994145762 , 1.000126065004428, 1.0263457672380667, 0.9982100819673314 , 1.000171493144087, 1.0307357727159434, 0.9979137265476951 , 1.000223866893561, 1.0351255884930834, 0.9976179336650787...
model Lorenz { const rho = 45.92 const beta = 4.0 const alpha = 16.0
state X, Y, Z
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 } }}
> T <- 1000> init_parameters <- list(X = 1, Y = 1, Z = 1)
> synthetic_dataset <- bi_generate_dataset(end_time=T, model=Lorenz, init=init_parameters)
> synthetic_data <- bi_read(synthetic_dataset)
> head(synthetic_data$X) time value1 0 1.000002 1 -12.653393 2 -13.725964 3 -13.573615 4 -11.268446 5 -15.36319
The end_time in LibBi is a real time, not a number of time steps. So it looks like your Haskell code simulates from time 0 to 0.001 with 1000 outputs, while your LibBi code simulates from time 0 to 1000 with 1000 outputs, which is bound to look different. To match your Haskell results in LibBi, you will want to set end_time to 0.001, and set the number of outputs to 1000 (I'm not sure what the exact argument for this is in RBi, with the command line interface it is --noutputs).
--
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.