# 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;
real sy = theta;
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 root;
real sx;
real sy;
real p;
sx = theta;
sy = theta;
p = theta;
root = integrate_1d(r_marginal_polar_bivariate_normal_pdf, 0, y, {sx, sy},
rep_array(0.0, 0), rep_array(0, 0), 1e-8)-p;
return root;
}
}
generated quantities {
real x_r;
int x_i;

vector y;
vector y_guess;
vector theta;
theta = 1;       // sx
theta = 1;       // sy
theta = .5;      // p
y_guess = 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);
}