Hi all,
I am learning to use Stan and attempting to build a model to find rings in a 2D image. These rings are circular in shape with a Gaussian profile in the radial direction, and an image can consist of a number of rings. We have written a model to fit to these images, by fitting to the central coordinates of the rings, the ring widths, the distance between the rings and the location of the first rings. We are enforcing that the distance between rings remains constant, so knowing the location of the first ring and the distance between rings is enough to specify the position of the rest of the rings.
There is a simple Python function below that we can use to generate the data:
def make_some_rings(x0, y0, xcoords, ycoords, radii, widths):
"""
Function to generate 15 rings on a 240, 80 2D image.
Args:
x0 (float): x-coordinate of centre of circle.
y0 (float): y-coordinate of centre of circle.
xcoords (1D-array): x-coordinates of 2D image pixels.
ycoords (1D-array): y-coordinates of 2D image pixels.
radii (list): radial value of each pixel in 2D image.
widths (list): standard deviation of each ring width in pixel space.
"""
# Generate radial position for every coordinate in the image
pixel_radii = np.sqrt((xcoords-x0)**2 + (ycoords-y0)**2)
image = np.zeros(xcoords.shape)
for i,dist in enumerate(radii):
# For each ring add a radial gaussian to each pixel
image += np.exp(-(pixel_radii-dist)**2/(2*widths[i]**2))
return image
We used 240x80 images and flattened them to a 1D list for analysis in our model. We adopt a statistical model we have previously found to work in other contexts. More specifically, we integrate over an uninformative prior on both the amplitude of the rings and the variance of the assumed Gaussian noise such that the likelihood is a Student-T.
Our Stan model is shown below:
data {
int<lower=1> N_d; // length of data vector
vector[N_d] pixels; // counts for all pixels
int<lower=1> N_r; // Number of rings
vector[N_d] x_pos; // x coordinate of every datapoint in 2D image
vector[N_d] y_pos; // y coordinate of every datapoint in 2D image
}
parameters {
real<lower=0,upper=1> delta_unit; // Distance between rings
real<lower=0,upper=1> r_1_unit; // Radius of first ring
real<lower=0,upper=1> x_0_unit; // x-coordinate of circle centre
real<lower=0,upper=1> y_0_unit; // y-coordinate of circle centre
vector<lower=0,upper=1>[N_r] sigma_r_unit; // Radial extent of each ring
}
transformed parameters {
real delta = (delta_unit*10) + 8;
real r_1 = (r_1_unit*30) + 5;
real x_0 = (x_0_unit*35) + 85;
real y_0 = (y_0_unit*20)+ 240;
vector[N_r] sigma_r = (sigma_r_unit*2) + 1;
}
model {
matrix[N_d, N_r] G_theta; //
vector[N_r] r;
real term_1;
real combined_2_3_4;
real term_2;
row_vector[N_r] term_3;
vector[N_r] term_4;
for (j in 1:N_r) {
r[j] = r_1 + (delta * (j - 1));
G_theta[:, j] = exp(-square(sqrt(square(x_pos - x_0) + square(y_pos - y_0)) - r[j]) ./ square(sigma_r[j]) * 2);
}
term_1 = 1 / sqrt(abs(determinant(G_theta' * G_theta)));
term_2 = dot_product(pixels', pixels);
term_3 = pixels'*G_theta * inverse(G_theta'* G_theta);
term_4 = G_theta' * pixels;
combined_2_3_4 = ((-N_d+N_r) / 2) * log(term_2 - term_3* term_4);
target += log(term_1) + combined_2_3_4;
// Prior
delta_unit ~ uniform(0, 1);
r_1_unit ~ uniform(0, 1);
x_0_unit ~ uniform(0, 1);
y_0_unit ~ uniform(0, 1);
sigma_r_unit ~ uniform(0, 1);
}
Issues we are encountering:
• Sampling and warm up are prohibitively slow, frequently complaining of over constrained parameters regarding G_theta. “ Exception: multiply: B[1] is nan, but must not be nan! (in ‘Stan_stage_1.stan’ at line 51)” Predominantly appears to happen in the warm up phase.
• Optimiser is inconsistent – is this due to there being many likelihood minima?
Any help with these issues would be greatly appreciated.
Thanks,
David