Let $n$ points $(x_i, y_i)\in \mathbb R^2$ satisfy $$ ax_i + bx_iy_i + cy_i = 1,\, i=1,\ldots, n $$ Each point is observed with an error $\boldsymbol \varepsilon _i = (\varepsilon _i^x, \varepsilon _i^y)\sim \mathcal N(\mathbf 0, \boldsymbol\Sigma _i)$, where $$ \boldsymbol\Sigma _i = \begin{bmatrix} \sigma _i^x & 0 \\ 0 & \sigma _i^y \end{bmatrix}^2 $$ and $(\boldsymbol \varepsilon _1,\ldots, \boldsymbol \varepsilon _n)$ are all independent. Given a dataset $$ \mathcal D = \begin{bmatrix} x_1 + \varepsilon _1^x & y_1 + \varepsilon _1^y & \sigma _1^x & \sigma _1^y \\ x_2 + \varepsilon _2^x & y_2 + \varepsilon _2^y & \sigma _2^x & \sigma _2^y \\ \vdots &\vdots &\vdots &\vdots \\ x_n + \varepsilon _n^x & y_n + \varepsilon _n^y & \sigma _n^x & \sigma _n^y \end{bmatrix}= \begin{bmatrix} \mathbf x^\varepsilon & \mathbf y^\varepsilon & \boldsymbol \sigma ^x & \boldsymbol \sigma ^y \end{bmatrix} $$ the goal is to estimate parameters $a, b, c$.

What have I already done? I expressed $y$ as a function of $x$ explicitly $$ y_i = \frac {-ax_i + 1}{bx_i + c} $$ and found a likelihood function: $$ -\ln p(\mathbf x^\varepsilon ,\mathbf y^\varepsilon \mid \mathbf x, \boldsymbol \sigma ^x, \boldsymbol \sigma ^y, a, b, c)\propto \sum _{i=1}^n\left[\frac {x^\varepsilon_i - x_i}{\sigma _i^x}\right]^2 + \left[\frac {y^\varepsilon_i - \frac {-ax_i + 1}{bx_i + c}}{\sigma _i^y}\right]^2 $$ which I want to minimize. For this purpose I wrote a simple R code:

generate_data <- function(n, a, b, c)
  x <- runif(n, 1, 15)
  y <- (-a*x + 1)/(b*x + c)

  sigma_x <- (rnorm(n, 0, 1))^2
  sigma_y <- (rnorm(n, 0, 2))^2

  data.frame(x=x, y=y, x_eps=x+rnorm(n, 0, sigma_x), y_eps=y+rnorm(n, 0, sigma_y), 
             sigma_x=sigma_x, sigma_y=sigma_y)

L <- function(theta, data)
  a <- theta[1]
  b <- theta[2]
  c <- theta[3]
  x <- theta[4:length(theta)]
  x_eps <- data$x_eps
  y_eps <- data$y_eps
  sigma_x <- data$sigma_x
  sigma_y <- data$sigma_y

  sum(((x_eps - x)/sigma_x)^2 + ((y_eps - (-a*x + 1)/(b*x + c))/sigma_y)^2)

data <- generate_data(50, 1000, 1, 10)
plot(x=data$x_eps, y=data$y_eps, pch="O")
points(x=data$x, y=data$y, col="red", pch=1)

optim(par=c(c(1, 1, 1), data$x_eps), fn=L, data=data, control=list(maxit=1e05))

However it does not converge to the optima. What can I do in this case?

  • $\begingroup$ The asymmetric treatment of $x$ and $y$ suggests this might not be a good approach. $\endgroup$ – whuber Nov 17 '15 at 13:35
  • $\begingroup$ @whuber Can you please explain your opinion? $\endgroup$ – nmerci Nov 17 '15 at 14:05
  • $\begingroup$ Consider the case where $a\approx 1$, $b\approx -1$, and $c\approx 1$, with values of $x$ close to $1$, with both positive and negative deviations. $y$, as a function of $x$, becomes singular within this region, making your approach extremely unstable. General mathematical meta-principles suggest that when a problem is formulated in a symmetric way, then its solution ought to respect that symmetry. That's certainly a good idea here. $\endgroup$ – whuber Nov 17 '15 at 14:16
  • $\begingroup$ That is true, but let us restrict to the case where $y$ (as a function of $x$) does not have discontinuities. The parameters in my R-example: $a=1000, b=1, c=10$ with $x\in (1, 15)$ will not "spoil" the example. Anyway, why the loss function stops converging so early? $\endgroup$ – nmerci Nov 17 '15 at 16:17
  • $\begingroup$ I wouldn't expect optim to be able to handle a problem of this size. You can minimize the loss much more effectively by finding the closest point on the curve to each data point (relative to the metric determined by the $\sigma^2$). That will reduce your problem to $3$ variables instead of $n+3$, which is an enormous improvement. $\endgroup$ – whuber Nov 17 '15 at 16:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.