CODEWITHSUNDEEP
roptimizationmaxstata
Francis Smart
Francis Smart

Reputation: 4055

Writing my own MLE command in R is causing issues

I am just really getting into trying to write MLE commands in R that function and look similar to native R functions. In this attempt I am trying to do a simple MLE with

y=b0 + x*b1 + u

and

u~N(0,sd=s0 + z*s1)

However, even such a simple command I am having difficulty coding. I have written a similar command in Stata in a handful of lines

Here is the code I have written so far in R. 

  normalreg <- function (beta, sigma=NULL, data, beta0=NULL, sigma0=NULL,
                         con1 = T, con2 = T) {

    # If a formula for sigma is not specified 
    #  assume it is the same as the formula for the beta.
    if (is.null(sigma)) sigma=beta

    # Grab the call expression
    mf <- match.call(expand.dots = FALSE)

    # Find the position of each argument
    m <- match(c("beta", "sigma", "data", "subset", "weights", "na.action", 
                 "offset"), names(mf), 0L)

    # Adjust names of mf
    mf <- mf[c(1L, m)]

    # Since I have two formulas I will call them both formula
    names(mf)[2:3] <- "formula"

    # Drop unused levels
    mf$drop.unused.levels <- TRUE

    # Divide mf into data1 and data2
    data1  <- data2 <- mf
     data1 <- mf[-3]
     data2 <- mf[-2]

    # Name the first elements model.frame which will be 
    data1[[1L]] <- data2[[1L]] <- as.name("model.frame")

    data1 <- as.matrix(eval(data1, parent.frame()))
    data2 <- as.matrix(eval(data2, parent.frame()))

    y     <- data1[,1]
    data1 <- data1[,-1]
     if (con1)  data1 <- cbind(data1,1)
    data2 <- unlist(data2[,-1])
      if (con2) data2 <- cbind(data2,1)

    data1 <- as.matrix(data1) # Ensure our data is read as matrix
    data2 <- as.matrix(data2) # Ensure our data is read as matrix

    if (!is.null(beta0)) if (length(beta0)!=ncol(data1))
      stop("Length of beta0 need equal the number of ind. data2iables in the first equation")

    if (!is.null(sigma0)) if (length(sigma0)!=ncol(data2)) 
      stop("Length of beta0 need equal the number of ind. data2iables in the second equation")

    # Set initial parameter estimates
    if (is.null(beta0))  beta0   <- rep(1, ncol(data1))
    if (is.null(sigma0)) sigma0 <- rep(1, ncol(data2))

    # Define the maximization function
    normMLE <- function(est=c(beta0,sigma0), data1=data1, data2=data2, y=y) {          
      data1est <- as.matrix(est[1:ncol(data1)], nrow=ncol(data1))
      data2est <- as.matrix(est[(ncol(data1)+1):(ncol(data1)+ncol(data2))],
                              nrow=ncol(data1))

      ps <-pnorm(y-data1%*%data1est, 
                       sd=data2%*%data2est)
      # Estimate a vector of log likelihoods based on coefficient estimates
      llk <- log(ps)
      -sum(llk) 
    }

    results <- optim(c(beta0,sigma0), normMLE, hessian=T,
                     data1=data1, data2=data2, y=y)

    results
  }


  x <-rnorm(10000)
  z<-x^2
  y <-x*2 + rnorm(10000, sd=2+z*2) + 10

  normalreg(y~x, y~z)

At this point the biggest issue is finding an optimization routine that does not fail when the some of the values return NA when the standard deviation goes negative. Any suggestions? Sorry for the huge amount of code.

Francis

Upvotes: 0

Views: 240

Answers (1)

Dason
Dason

Reputation: 61933

I include a check to see if any of the standard deviations are less than or equal to 0 and return a likelihood of 0 if that is the case. Seems to work for me. You can figure out the details of wrapping it into your function.

#y=b0 + x*b1 + u
#u~N(0,sd=s0 + z*s1)

ll <- function(par, x, z, y){
    b0 <- par[1]
    b1 <- par[2]
    s0 <- par[3]
    s1 <- par[4]
    sds <- s0 + z*s1
    if(any(sds <= 0)){
        return(log(0))
    }

    preds <- b0 + x*b1

    sum(dnorm(y, preds, sds, log = TRUE))
}

n <- 100
b0 <- 10
b1 <- 2
s0 <- 2
s1 <- 2
x <- rnorm(n)
z <- x^2
y <- b0 + b1*x + rnorm(n, sd = s0 + s1*z)

optim(c(1,1,1,1), ll, x=x, z=z,y=y, control = list(fnscale = -1))

With that said it probably wouldn't be a bad idea to parameterize the standard deviation in such a way that it is impossible to go negative...

Upvotes: 2

Related Questions