How to specify operations on parameters? (Damped Oscillator model)

[rohitbatra29]

Hi,

I'm trying to fit a damped oscillator model in ctsem and facing issues:

$$d^{2}y_{t} = \left( \omega \cdot (y_{t} - \mu) + \zeta \cdot \frac{dy_{t}}{dt} \right) dt^{2} + \sigma \cdot dW_{t}$$

I transformed it into the following two equations for level and velocity:

$$ \begin{align} dx_{t} &= p_{t} \\ dp_{t} &= \omega \cdot (x_{t} - \mu) + \zeta \cdot p + \sigma \cdot dW_{t} \end{align} $$

Following is the code I used:

univariate_oscillator <- ctModel(
  type = "stanct",
  n.latent = 2, 
  n.manifest = 1,
  manifestNames = c("ibi"),
  latentNames = c("X", "P"),
  time = "time",
  id = "id",
  LAMBDA = matrix(c(1, 0), nrow = 1),
  DRIFT = matrix(c(0, 1, 
                   "omegax", "zetax"), nrow = 2, byrow = TRUE),
  CINT = matrix(c(0, 
                  "-omegax * mu_x"), nrow = 2, byrow = TRUE), # it works if I just estimate one parameter instead of -omegax*mu_x
  DIFFUSION = matrix(c(0, 0, 
                       0, "sigmax"), nrow = 2, byrow = TRUE),
  MANIFESTMEANS = diag(0, 1),
  MANIFESTVAR = diag(0, 1),
  T0MEANS = matrix(c(mu_x0, 0), nrow = 2, byrow = TRUE),
  T0VAR = diag(1, 2),
  PARS = c("mu_x")
)

univariate_oscillator$pars$indvarying <- FALSE


ctModelLatex(univariate_oscillator)

fit_mod <- ctStanFit(datalong = p1_data, ctstanmodel = univariate_oscillator)

I'm facing issues with the operations on parameters and was wondering if you had any suggestions on how to code this model? The error I get is:

> fit_mod <- ctStanFit(datalong = p1_data, ctstanmodel = univariate_oscillator)
Maximum likelihood estimation requested
Setting DRIFT diagonals of 0 to approximate 0 of -1e-06
Compiling model...
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Info: integer division implicitly rounds to integer. Found int division: d * d - d / 2
 Positive values rounded down, negative values rounded up or down in platform-dependent way.
Variable "NA" does not exist.
 error in 'model3ad63ff9f2f_ctsem' at line 202, column 79
  -------------------------------------------------
   200: 
   201: if(multiplier != 1.0) param *=multiplier;
   202: if(transform < 49 && offset != 0.0) param+=offset;if(transform==-11) param = NA;
                                                                                      ^
   203:     return param;
  -------------------------------------------------

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'ctsem' due to the above error.

My impression is that the way I'm specifying the parameter multiplication is not being read correctly for stan's syntax. But I also tried running the following example from the manual (which looked to have an operation on parameters) but even this model did not run for me:

sunspots<-sunspot.year
sunspots<-sunspots[50: (length(sunspots) - (1988-1924))]
id <- 1
time <- 1749:1924
datalong <- data.frame(id, time, sunspots)

m <- ctModel(type='stanct',
manifestNames='sunspots',
latentNames=c('ss_level', 'ss_velocity'),
LAMBDA=c( 1, 'ma1|log(1+exp(param))'),
DRIFT=c(0, 1,
'-log1p_exp(freqintercept + freqbystate * ss_level)','a22'),
MANIFESTMEANS=c('m1|param * 10 + 44'),
MANIFESTVAR=diag(0,1), #As per original spec
CINT=0,
DIFFUSION=c(0, 0,
0, "diffusion"),
PARS=c('freqintercept', 'freqbylevel'))

m$pars$indvarying <- FALSE

ssfitnl <- ctStanFit(datalong, m, optimcontrol=list(stochastic=TRUE))

[cdriveraus]

I don't think you need a nonlinear model for a damped linear oscillator -- the manual ( https://cran.r-project.org/web/packages/ctsem/vignettes/hierarchicalmanual.pdf ) shows a linear version with the sunspots data. To proceed with the nonlinear form, you need to make sure all your 'complex' parameters are specified as additional parameters in the PARS matrix -- you've only specified mu_x , but you also need to specify omega_x as an extra parameter. The example you found that is not working doesn't exist in the current ctsem code, I guess there was a temporary bug at some point -- to fix it you just need to make sure the PARS are properly specified, change 'freqbystate' in the DRIFT matrix to 'freqbylevel' which we have an additional parameter for.

[rohitbatra29]

Great, that did the trick! I assumed we wouldn't need to to declare omega_x since it's already declared elsewhere. Thank you!

[cdriveraus]

Great :)

…

On Sat, Jan 18, 2025, 00:25 Rohit Batra ***@***.***> wrote: Great, that did the trick! I assumed we wouldn't need to to declare omega_x since it's already declared elsewhere. Thank you! — Reply to this email directly, view it on GitHub <#62 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA7LCBI775GJ3Z4TNOAM4O32LGGPPAVCNFSM6AAAAABVKQH5DSVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCOBXGIZTMOA> . You are receiving this because you commented.Message ID: ***@***.***>