A term in my log-likelihood can only be numerically approximated. To approximate this term, Brent’s root-finding method (which corresponds to function `uniroot`

in base R) needs to be used. Since the Stan Math library does not have the counterpart of `uniroot`

, I am planning to use an external C++ function that is not in the Stan Math library to perform the root-finding step.

My questions are :

a. Would numerical approximation of the term using a root-finding method prevent the autodifferentiation of Stan from working?

b. How should I import the external C++ function that is not in the Stan Math library? Do I simply follow the steps in interfacing with external C++ code (https://cran.r-project.org/web/packages/rstan/vignettes/external.html)?

Thank you!

1 Like

Numerical methods are fine, but the root is an implicit function. So, you have to use the implicit function theorem to get the correct derivatives with respect to the unknown parameters. I would start with `algebraic_solver`

in Stan before trying to use `toms_748`

from Boost, but if you were to use `toms_748`

, it would be like in that vignette.

3 Likes

Thank you for your prompt response. I used `algebra_solver`

as you suggested to define the following function to compute one term in my log likelihood function. Unfortunately, I ran into a parser error that says PARSER EXPECTED: “,”. I understand it is caused by not passing function arguments x_r and x_i, so I attempted to add line a and b (see the commented code below) and pass x_r and x_i to the function (see line c). Then the parser says “fourth argument to algebra_solver must be data only (cannot reference parameters).have to be data (not parameters)." In my case, data is not necessary to solve the equation. I am not sure how to handle this?

Below is the relevant code:

```
real cal_some_term(vector gam, vector lam) {
real res;
int p = num_elements(gam);
vector[p] lam2 = lam;
// line a: real x_r[0];
// line b: int x_i[0];
// some other code
vector[2*p] theta;
vector[1] guess;
theta = append_row(gam, lam2);
guess[1] = (low + up) / 2.0;
tau = algebra_solver(system, guess, theta);
// line c: tau = algebra_solver(system, guess, theta, x_r, x_i);
// code to calculate res as a function of tau
return res;
}
```

Thank you!

You need to specify `x_r`

and `x_i`

even if they are of size zero and not used by your `system`

function. Nevertheless, the `system`

function must accept those two arguments.

1 Like

So is it sufficient to put `x_r`

and `x_i`

as dummy data in the data section? e.g.

```
data {
int<lower=0> N;
int<lower=0> K;
vector[K] y[N];
real x_r[0];
int x_i[0];
}
```

and pass any value to `x_r`

and `x_i`

.

Almost – you just need to move the empty arrays into the `transformed data`

block,

```
data {
...
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
...
}
```

2 Likes