# Algebra_solver with integrate_1d containing modified_bessel_first_kind

I created a function for the PDF of an uncorrelated bivariate normal distribution as a function of radius (equation 17 from A Toolbox for the Radial and Angular Marginalization of Bivariate Normal Distributions (arxiv.org)). My goal is to integrate that, and solve to return the quantile associated with a given cumulative probability. For example, in the code below I want to find the radius r that contains p = 50\% of the data given an uncorrelated bivariate normal distribution centered at the origin with standard deviations \sigma_x = 1, and \sigma_y = 1.

I have tested this and can integrate with no problem, but trying to solve for p is a problem that seems to be due to modified_bessel_first_kind. If I comment out the line with modified_bessel_first_kind it works, but otherwise I get the error: Exception: gradient_of_f: The gradient of f is nan for parameter 0. Is there any way that I might be able to work around this?

functions{
real r_marginal_polar_bivariate_normal_pdf(
real x, real xc, real[] theta, real[] x_r, int[] x_i) {

// define my parameters
real sx = theta[1];
real sy = theta[2];
real prob;

prob = (x/(sx*sy)) *
exp(-(sx^2 + sy^2)/((2*sx*sy)^2) * x^2) *
modified_bessel_first_kind(0, -(sx^2 - sy^2)/((2*sx*sy)^2) * x^2 );
return prob;
}

vector system(vector y, vector theta, real[] x_r, int[] x_i){
vector[1] root;
real sx;
real sy;
real p;
sx = theta[1];
sy = theta[2];
p = theta[3];
root[1] = integrate_1d(r_marginal_polar_bivariate_normal_pdf, 0, y[1], {sx, sy},
rep_array(0.0, 0), rep_array(0, 0), 1e-8)-p;
return root;
}
}
generated quantities {
real x_r[0];
int x_i[0];

vector[1] y;
vector[1] y_guess;
vector[3] theta;
theta[1] = 1;       // sx
theta[2] = 1;       // sy
theta[3] = .5;      // p
y_guess[1] = 1.17;

// find the quantile of r corresponding with p
y = algebra_solver(system, y_guess, theta, x_r, x_i, 1e-8, 1e-5, 1e4);
}