# Convergence issues for fitting a normal linear measurement error model with optim in R

I've been trying to run a simulation on this normal theory linear measurement error model as provided by Boos and Stefanski Essential Statistical Inference pg. 114 but am having computational errors in the estimation. The model is as follows

$$Y_i\sim N(\alpha+\beta U_i, \sigma_\epsilon^2)$$

$$X_i \sim N(U_i, \sigma^2)$$

We observe $$\bf{Y}$$ and $$\bf{X}$$ for $$i=1,...,n$$, and the parameters of interest are $$\theta=(\alpha, \beta, \sigma_\epsilon^2)$$. The data are independent and $$\sigma^2$$ is known. However we have nuisance parameters $$U_1,...,U_n$$ that grow at the same rate as the sample size $$n$$. The full likelihood $$L(\theta, \textbf{U}|\textbf{Y,X})$$ diverges and we are unable to conduct standard MLE in this scenario.

The strategy is to use conditional likelihood. The book provides a sufficient "statistic" $$T_i=Y_i\beta / \sigma_\epsilon^2+X_i / \sigma^2$$ for $$U_i$$. Then the distribution $$[\textbf{Y}| \theta, \textbf{U},\textbf{T}]=[\textbf{Y}|\theta, \textbf{T}]$$ is independent of the nuisance parameters and we can use this conditional likelihood to estimate the parameters $$\theta$$. The contribution to the conditional likelihood of an individual observation is proportional to the distribution $$Y|T$$

$$Y_i|T_i=t\sim N\left((\alpha/\sigma^2+\beta t)\frac{\sigma_\epsilon^2\sigma^2}{\sigma^2\beta^2+\sigma^2_\epsilon}, \frac{\sigma^4_\epsilon}{\sigma^2\beta^2+\sigma^2_\epsilon}\right)$$

The code I'm using to implement this strategy is below. However the estimation is way off and seems very unstable. I haven't been able to figure out where things could be going wrong?

set.seed(1)

# knowns
n     <- 100
sigma <- 10

# unknowns
alpha   <- 2
beta    <- 4
sigma_e <- 2
U       <- runif(n, 0, 1)

# Observed data
X  <- rnorm(n, U, sigma)
Y  <- rnorm(n, alpha + beta*U, sigma_e)

# f(Y|T=t)
log_conditional_pdf <- function(y, t, alpha, beta, sigma_e, sigma){
mu_cond <- (alpha / sigma^2 + beta*t) * sigma_e^2 * sigma^2 / (sigma^2 * beta^2 + sigma_e^2)
sigma_cond <- sqrt( sigma_e^4 / (sigma^2 * beta^2 + sigma_e^2) )
dnorm(y, mu_cond, sigma_cond, log = TRUE)
}

loss <- function(theta, Y, X, sigma){
alpha   <- theta[1]
beta    <- theta[2]
sigma_e <- exp(theta[3])

# The sufficient "statistic"
T_i <- Y * beta / sigma_e^2 + X / sigma^2

loglik <- log_conditional_pdf(Y, T_i, alpha, beta, sigma_e, sigma)

return(-sum( loglik ))
}

########## Results #####################################

optim(c(1, 1, 2), loss, sigma = sigma, Y = Y, X = X)
# $par # [1] 3.949812 10.181010 -15.223716 #$value
# [1] -3357.989

# $counts # function gradient # 410 NA #$convergence
# [1] 0

optim(c(alpha, beta, log(sigma_e)), loss, sigma = sigma, Y = Y, X = X)

# $par # [1] 4.301806 12.868783 -15.183607 #$value
# [1] -3370.137

# $counts # function gradient # 331 NA #$convergence
# [1] 10

$$`$$