search

traitecoevo / hmde

Implementation of hierarchical Bayesian longitudinal models to estimate differential equation parameters.
Other
2 stars 3 forks source link

Changes to numerical integration #32

Open Tess-LaCoil opened 2 months ago

Tess-LaCoil commented 2 months ago

In order to avoid pathologies associated with particular numerical methods:

  1. Use analytic solutions for von Bertalanffy, power-law (log transformed)
  2. Use the inbuilt stiff solver for Canham
  3. Retain RK4 for linear model for demonstration purpose
  4. Construct a function that tests stability of numerical methods for posterior parameter estimates
Tess-LaCoil commented 2 months ago

@dfalster can you please have a look at the implementation of models in the intergration-update branch and check if you believe the maths.

So far I have done analytic solutions for von Bertalanffy based on $$Y(t) = Y{max} + (Y(0) - Y{max}) \exp(-\beta t)$$ and a log-transformation of the power law with $$\log(Y(t)) = \frac{\log(\beta_0)}{\beta_1} + \Bigg(\log(Y(0)) - \frac{\log(\beta_0)}{\beta_1}\Bigg) \exp(-\beta t)$$ where $\beta_0$ is the coeff parameter, $\beta_1$ the power.

For the last estimate of growth when there's no following measurement to compare I have it take the same time gap as the previous observation period and project forward, then use that difference.

Tess-LaCoil commented 2 months ago

von Bertalanffy model is fully updated. Now working on power law.

Tess-LaCoil commented 2 months ago

Numerical solution with step size 0.1 NumericSolution

Analytic solution AnalyticSolution In wrangling the power law model I have come across what I think is an error in the maths. I'm going to ask a friend to double check my solution but the RK4 model reliably gives the correct sizes over time and parameters while the analytic solution does not and I have been unable to determine why.

The current models are on the integration-update branch, implemented as separate stan files in the inst/stan/ directory. @dfalster can you please look at the power_single_ind_analytic.stan file and see if you believe the maths? To run it, save a duplicate under the name power_single_ind.stan and then use the usual workflow.

Tess-LaCoil commented 2 months ago

After some additional wrangling and checking of the solution, what breaks the numerical methods for vB doesn't appear in the power law model because the coefficient on Y is not negative. It may be preferable to sidestep the problems in the analytic solution by just using the numerical method instead as it reliably gives good estimates.

Tess-LaCoil commented 2 months ago

image Power law parameter estimates are robus to changes in the step size.