Bayesian lighthouse location estimation

Question

I am trying to learn Stan in R and as a fun challenge I am trying to estimate the location of a lighthouse based on the observed flashes. But the models I tried do not converge (Rhat != 1) or have estimated parameters with a large spread.

The observed data are flashes from a lighthouse 100 meters away the (straight) coast line. The angle is uniformly distributed, but the observed flashes along the coast line are heavy tailed.

n_flashes <- 50
loc <- c(0, 100)

angles <- runif(n_flashes, -pi/2, pi/2)
angles_x <- loc[2] * tan(angles)
flashes <- loc[1] + angles_x

I want to estimate the lighthouse location by given the actual model generation process to stan. This is my stan model:

data {
  int<lower=0> N;
  real flashes[N];
}

parameters {
  real x_loc;
  real<lower=0> y_loc;
  real<lower=-pi()/2, upper=pi()/2> angle[N];
}

model {
  x_loc ~ normal(0, 10);
  y_loc ~ normal(100, 10);

  for (i in 1:N) {
    flashes[i] ~ cauchy(x_loc + tan(angle[i]) * y_loc, 1);
  }
}

(Note: I model the flash observation with a 1 meter standard deviation, because Stan requires me to give a distribution for flashes[i]).

Then I call the model from R:

n_flashes <- 50
loc <- c(0, 100)

angles <- runif(n_flashes, -pi/2, pi/2)
angles_x <- loc[2] * tan(angles)
flashes <- loc[1] + angles_x

stan_input <- list()
stan_input$flashes <- flashes
stan_input$N <- n_flashes
fit <- stan("lighthouse.stan", data = stan_input)

But this model has large Rhat's. How can I improve this model to prevent this? How can I model the angles in stan? How would you model the location of a lighthouse based on the flashes observed along the coast line?

Answer 1

This is a famous problem known as Gull's Lighthouse, from an example by Gull in 1988. It has deep implications when taken one additional step in both the social sciences and in physics. You actually have enough information to solve this problem through acceptance-rejection testing, but if you want to use MCMC feel free.

Let's look at your problem another way, knowing that the Lighthouse is 100 m from shore is actually a lot of information. That piece is a big deal, otherwise you would have to solve both distance and location.

As a side note, the Cauchy distribution is 636 semi-interquartile ranges wide at 99.99% HDR so if you were looking for narrow intervals, you can almost forget it. If you use a normal approximation, however, it will never converge as this problem violates the central limit theorem and if you estimate $\sigma$ as standard deviation then if $n$ is sample size then $s^2\to\infty$ as $n\to\infty$. If you did that, your estimates would get worse and worse as you added date points.

You need to note that you are observing the set $\{x_i,i\in{1\dots{n}}\}$ and there is some point $\mu$ where the lighthouse would be perpendicular to the beach. Note that it is not the positioning of the $x_i$ which are uniformly distributed, but rather the angle, $\theta$. We did not observe the angle, only the set of points.

The relationship of the points to the angle is $$\frac{x_i-\mu}{100}=\tan\theta.$$ In terms of the angle, this can be expressed as $$\theta=\tan^{-1}\left(\frac{x_i-\mu}{100}\right)$$

Since $$p(x_i|\mu,100)=p(\theta|\mu,100)|\frac{\mathrm{d}\theta}{\mathrm{d}x_i}|$$ and $\theta$ is uniform over the region $(-\pi/2,\pi/2)$ we know that $$p(x_i|\mu,100)=\frac{1}{\pi}|\frac{\mathrm{d}\theta}{\mathrm{d}x_i}|,$$ which is $$p(x_i|\mu,100)=\frac{1}{\pi}\frac{100}{100^2+(x_i-\mu)^2}.$$

So, by Bayes theorem, if the shore point perpendicular to the lighthouse is $(\mu,0)$, then the lighthouse is located at $(\mu,100)$. The estimator of this point, since you have no prior information on $\mu$ is $$p(\mu|x_1\dots{x_n},100)\propto\prod_{i=1}^n\frac{1}{\pi}\frac{100}{100^2+(x_i-\mu)^2},\forall\mu\in\Re.$$

If you want to make the problem more interesting, assume you do not know that the location is one hundred meters from shore. Your problem them becomes the two dimensional problem more suited to stan. Your problem then becomes $$p(\mu,\gamma|x_1\dots{x_n})\propto\prod_{i=1}^n\frac{1}{\pi}\frac{\gamma}{\gamma^2+(x_i-\mu)^2},\forall\mu\in\Re;\gamma\in\Re^{++}.$$

As a side note, this problem is related to the also famous Witch of Agnesi, studied by Maria Agnesi and published in 1748 and the more basic problem studied by Fermat and Grandi in 1703. She had included it in the very first true calculus textbook. It provided students an education from algebra through integral and differential equations.

Answer 2

The flashes' positions are modeled precisely by the Cauchy, as @user25459 said, so it's no coincidence that Stan is able to estimate the true values using this model (I use $\alpha$ for $\mu$, and $\beta$ for $\gamma$):

model <- '
  data {
    int<lower=0> N;
    real x_[N];
  }

  parameters {
    real alpha;
    real<lower=0> beta;
  }

  model {
    alpha ~ uniform(0, 20);
    beta  ~ uniform(0, 50);

    for (k in 1:N) {
      x_[k] ~ cauchy(alpha, beta);
    }    
  }
'

Let's try this in Stan:

alpha <- 10 # unknown true values
beta  <- 30 
##################
set.seed(123)
N       <- 100
theta_k <- runif(N,-pi/2,pi/2)
x_k     <- beta * tan(theta_k) + alpha

stan_input <- list(x_=x_k, N=N)

fit  <- stan(model_code=model, data=stan_input, iter=1000, verbose=FALSE)
fit2 <- stan(fit=fit, data=stan_input, iter=5000, warmup=2000, verbose=FALSE)
print(fit2, pars=c("alpha", "beta"))

Get us:

Inference for Stan model: 36eef3c3637fb3a9564529926f8463fe.
4 chains, each with iter=5000; warmup=2000; thin=1; 
post-warmup draws per chain=3000, total post-warmup draws=12000.

       mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
alpha  8.96    0.05 3.91  1.61  6.13  8.95 11.68 16.80  7079    1
beta  30.29    0.05 4.18 22.84 27.45 29.98 32.88 39.18  8470    1

Samples were drawn using NUTS(diag_e) at Wed Jan 11 18:28:57 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Let's plot one of the samples:

la <- extract(fit2)
alpha_hats <- as.numeric(la$alpha)  # vertical line is true value

About the angles, notice that

$$\theta_k = \arctan \frac{x_k - \alpha}{\beta}$$

Each $\theta_k$ is a function of two random vars, so is itself a random var.

To find the empirical distribution of angle $\theta_k$ we use the random samples of $\alpha$ and $\beta$.

theta_1 <- atan((x_k[1] - alpha_hats)/beta_hats)
hist(theta_1, breaks=100, xlab=bquote(theta[1]), ylab="", main="Distribution for 1st angle", yaxt='n')