r - Rでシミュレートされた可能性を最大化します - 私が得ていない明らかなものはありますか?

Question

私は、周波数をシミュレートし、それらを確率として使用することによって、離散選択 (Lerman and Manski (1981)) でシミュレートされた可能性を最大化しようとしています (直接計算することはできません)。ただし、R は最適な値を見つけることができません (最大化は常に開始値を生成します)。最小限の例として、非常に単純なプロビット推定のコードを次に示します。

### simulate data
set.seed(5849)
N <- 2000
b.cons  <- 8 
b.x         <- 10
x <- cbind(rep(1, N), runif(N)) #"observed variables"
e <- rnorm(N) # "unobserved error"
k <- runif(N)*10+7 # threshold: something random, but high enough to guarantee some variation in i
t <- x%*%c(b.cons, b.x)+e
i <- 1*(k>t) #participation dummy


### likelihood function
R   <- 1000 # number of draws
err <- matrix(rnorm(R*N), N, R) # draw error terms (outside of likelihood function to speed up estimation)

# estimate b.i, sig.i
probit.sim <- function(params, I, K, X) {
    part  =matrix(NA, N, R)
    T       = X%*%params%*%rep(1, R) + err
    for (i in 1:R) part[,i] = K>T[,i]

    pr.i1       = rowSums(part)/R
    pr.i1[pr.i1==0] <- 0.001 
    pr.i1[pr.i1==1] <- 0.999
    pr.i0       = 1-pr.i1

    llik = t(I)%*%log(pr.i1) + t(1-I)%*%log(pr.i0)
    -llik
}

### maximize likelihood
optim(c(1,1), probit.sim, I = i, K = k, X = x)

確率が十分に滑らかでないためですか？超スムーズでないものを最大化する方法はありますか? グラフでは、最大値はまだかなりはっきりしているように見えます...または、何か他のものを完全に見逃していますか?

初心者なので色々とアドバイスいただけると助かります！

(また、そのようなシミュレートされた最大尤度関数をプログラムする方法の詳細に実際に入る参考文献-私が見たほとんどの参考文献は、それについて非常に理論的なままです)