R에 내 MLE 명령을 쓰면 문제가 발생합니다.

Question

저는 R에 MLE 명령을 쓰려고 시도하고 있으며 R 함수와 비슷하게 보입니다. 이러한 시도에서, I는

Y로 간단한 MLE를 수행하는 것을 시도하고있다 = 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)) 

는`) (

이