Noisy samples from a disk of unknown radius
Suppose that we obtain “noisy” measurements of points sampled uniformly from a circular area of unknown radius. Our goal will be to infer the radius of the circular area, \(R\).
Suppose further that we know in advance:
- The position of the center of the circular area; in this case it is located at the origin, \((0, 0)\).
- Points that are uniformly sampled from a circular area with radius \(R\) induce a uniform distribution on the squared radius \(R^2\). That is, at any given angle \(\omega\), the points lying along the line from \((0, 0)\) to \((R \cos\omega, R \sin \omega)\) have uniformly distributed squared radial magnitudes.
- The noise terms act additively on the squared radii; this is a bit contrived, but it will make the modeling easier, and it implies that at any given angle \(\omega\), the noise perturbations will act radially outward or inward.
- The form of the noise distribution; in this case \(N(0, \sigma^2)\), with \(\sigma^2\) unknown.
The second and third items on the list above tells us that if we observe data drawn uniformly from the surface of a circle with radius \(R\), each point’s squared distance from the center of the circle will be uniformly distributed on \((0, R)\). Given these facts, let’s think about how we might do inference on \(R\).
0. Visualize the data generating process
It is often helpful to have a picture of what’s going on when we describe a generative model. Here we can simulate and plot our data, overlaying the true radius \(R\) on the sampled points \(r_1, \dots, r_n\):
#
n <- 1000
true_Rsquared <- 5
true_sigma <- 1.25
#
u <- runif(n, 0, true_Rsquared)
r <- u + rnorm(n, sd = true_sigma)
theta <- runif(n, 0, 2*pi)
#
ggplot2::ggplot() +
geom_point(data = data.frame(x = sign(r)*sqrt(abs(r))*cos(theta), y = sign(r)*sqrt(abs(r))*sin(theta)),
aes(x = x, y = y), shape = 1) +
geom_path(data = data.frame(R = true_Rsquared) %>%
plyr::ddply(.(R), function(d){
data.frame(x = sqrt(d$R)*cos(seq(0, 2*pi, length.out = 100)),
y = sqrt(d$R)*sin(seq(0, 2*pi, length.out = 100)))
}),
aes(x = x, y = y), alpha = 1, colour = "red") +
coord_fixed()
1. Write down the rest of the generative model
Let \(r_i\) denote the noisy squared radius of observed data point \(i\). Given the noiseless squared radial positions of the points, we assume that the observations \(r_i\) are normally distributed, each with mean \(u_i\) and variance \(\sigma^2\). A model for the data generating process might look something like this: \[ \begin{aligned} r_i | \mu_i, \sigma^2 &\sim N(u_i, \sigma^2) \\ u_i | R^2 &\overset{iid}\sim U(0, R^2), \\ \end{aligned} \] with priors: \[ \begin{aligned} R^2 | m,k &\sim Pa(m, k) \\ 1 / \sigma^2 | \alpha, \beta &\sim Ga(\alpha, \beta) \end{aligned} \]
Here we have used the conjugate prior, the Pareto distribution, for the uniform likelihood. We have also used the conjugate Gamma prior for the precision parameter.
2. Write down the data likelihood
Write \(r = (r_1, \ldots, r_n)\). The generative model implies that the likelihood takes the form
\[ \begin{aligned} L(r ; u, \sigma^2) &= \left(\frac{1}{2 \pi \sigma^2}\right)^{n/2} \textrm{exp} \left\{-\frac{1}{2 \sigma^2}\sum_{i=1}^n (r_i - u_i)^2\right\} \end{aligned} \]
3. Write down (and factorize) the joint distribution for all parameters
Which parameters are independent of one another? How does this inform the way we can decompose the joint distribution of the parameters and the data? What does this decomposition imply for a Gibbs sampling algorithm that produces draws from the posterior distribution of \(R^2\)?
\[ \begin{aligned} p(r, u, 1/\sigma^2, R^2 | m, k, \alpha, \beta) &= p(r | u, 1/\sigma^2, R^2, m, k) \cdot p(u | R^2) \cdot p(R^2 | m, k) \cdot p(1/\sigma^2 | \alpha, \beta) \end{aligned} \]
4. Write down (if possible) full conditional distributions
In this case, we can explicitly write out the full conditional distributions for the model parameters \(u, \theta\) and \(1/\sigma^2\). What are they? Make sure to keep track of indicator functions that may truncate the support of a density function.
\[ u_i | r, \sigma^2, R^2 \sim \text{Truncated Normal}\left(r_i, \sigma^2, 0, R^2 \right) \]
\[ 1 / \sigma^2 | r, u, R^2 \sim Ga\left(n/2 + \alpha, \frac{1}{2}\sum_{i=1}^n (r_i - u_i)^2 + \beta \right) \]
\[ R^2 | u, r, \sigma^2 \sim Pa\left(\max(u_1, \dots, u_n, m), k + n \right) \]
If you are seeing the truncated normal distribution for the first time, take a few minutes to read https://en.wikipedia.org/wiki/Truncated_normal_distribution.
