Would numerical approximation of log-likelihood prevent the autodifferentiation of Stan from working?

Interfaces RStan
1

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
2

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
3

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!

4

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
5

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.

6

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