CODEWITHSUNDEEP
r
user3443183
user3443183

Reputation: 115

Optimization has not converged in R

I have a function which estimates the two moments of a truncated Gaussian distribution by using the maximum likelihood method. In my database, I do not always obtain results which converge. So I put a warning to see how many of them they are.

Then I tried to erase them, by using the tryCatch function but it seems that I don't know how to use it properly.

Does someone have an idea to do this ?

Here is the code :

df=as.data.frame(matrix(0,1000,2))
df$V1=runif(nrow(df),0,1)
df$V2=sample(c(1:10),nrow(df), replace=TRUE)

library(truncnorm)
f <- function(d) # d is a vector
{
tryCatch(
{
f2 <- function(x) -sum(log(dtruncnorm(d, a=0, b=1, mean = x[1], sd = x[2])))
  res <- optim(par=c(mean(d),sd(d)),fn=f2)
  if(res$convergence!=0) warning("Optimization has not converged")
  return(list(res1=res$par[1],res2=res$par[2]^2))
  },
error = function(error) {return()} 
  )
}

res1=tapply(df$V1, df$V2, function(x) f(x)$res1)

So the idea is to "clean" res1 by deleting values which have not converged.

Best regards.

Upvotes: 0

Views: 1364

Answers (1)

Roman Luštrik
Roman Luštrik

Reputation: 70643

Here is a small example of how tryCatch works. You need to run your expression within tryCatch and look at the result. If the expression encounters an error (or a warning), you can adjust the workflow. The point is, that even though a warning or error are produced, life goes on and the program doesn't stop.

my.list <- vector("list", 15)

set.seed(357)
for (i in 1:length(my.list)) {
  result <- tryCatch(
    {
      if (runif(1) > 0.5) {
        stop("More than 0.5")
      } else {
        42
      }
      },
    error = function(e) e
    )

  # This bit catches the result and manipulates it further (output NAs or whathaveyou)
  if (!is.numeric(result)) {
    my.list[[i]] <- "Failed to converge, coerce to NA"
  } else {
    my.list[[i]] <- result
  }
}

my.list

[[1]]
[1] 42

[[2]]
[1] 42

[[3]]
[1] 42

[[4]]
[1] 42

[[5]]
[1] "Failed to converge, coerce to NA"

[[6]]
[1] 42

[[7]]
[1] "Failed to converge, coerce to NA"

[[8]]
[1] "Failed to converge, coerce to NA"

[[9]]
[1] "Failed to converge, coerce to NA"

[[10]]
[1] "Failed to converge, coerce to NA"

[[11]]
[1] "Failed to converge, coerce to NA"

[[12]]
[1] "Failed to converge, coerce to NA"

[[13]]
[1] 42

[[14]]
[1] 42

[[15]]
[1] "Failed to converge, coerce to NA"

If the expression doesn't evaluate (it encounters stop) the error message is returned. You can substitute this for your own function that outputs a data.frame with NAs or has some other desired behavior.

# custom function that produces custom output
myFun <- function(x) {
  NA
}

set.seed(357)
for (i in 1:length(my.list)) {
  result <- tryCatch({
    if (runif(1) > 0.5) {
      stop("More than 0.5")
    } else {
      42
    }
  },
  error = myFun
  )

  my.list[[i]] <- result

}

my.list

[[1]]
[1] 42

[[2]]
[1] 42

[[3]]
[1] 42

[[4]]
[1] 42

[[5]]
[1] NA

[[6]]
[1] 42

[[7]]
[1] NA

[[8]]
[1] NA

[[9]]
[1] NA

[[10]]
[1] NA

[[11]]
[1] NA

[[12]]
[1] NA

[[13]]
[1] 42

[[14]]
[1] 42

[[15]]
[1] NA

Upvotes: 1

Related Questions