在R中寫入我自己的MLE命令導致問題

Question

我只是真的嘗試在R中編寫MLE命令，它的功能與外部R函數類似。在這種嘗試我試圖做一個簡單的MLE與

Y = B0 + X * B1 + U

和

U〜N（0，SD = S0 + Z * S1）

但是，即使是這樣一個簡單的命令，我也很難編碼。我在Stata in a handful of lines

這裏寫了一個類似的命令是我在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) 

至今寫在這一點上，最大的問題是找到一個優化程序不失敗時，一些代碼當標準偏差變爲負值時，這些值返回NA。有什麼建議麼？對不起，代碼量很大。

Francis

Answer 1

我包括檢查以查看是否有任何的標準偏差均小於或等於0並返回似0程度，如果是這樣的話。似乎爲我工作。你可以弄清楚將它包裝到你的函數中的細節。

#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)) 

隨着中說，它可能不會是一個壞主意參數以這樣的方式，這是不可能去負的標準差......