2

我是一个 R 菜鸟，这可能反映在不那么密集的代码中 - 所以请忍耐。我正在尝试使用最大值估计二元正态分布的系数。似然估计。我在调用 OPTIM 函数时收到与 Hessian 相关的错误。我已经尝试了很多调试，但似乎无法摆脱错误。非常感谢您对如何解决这个问题的任何见解。

我使用的数据是 {y1,y2,x1,x2} 其中 y1,y2 是二进制变量。我用来模拟数据的代码如下：

x1=rnorm(1000)*2+3
x2=rnorm(1000)-0.5

mu=c(0,0)
sigma=array(c(1,0.5,0.5,1),c(2,2)) # correlation matrix
e=mvrnorm(n = 1000, mu, sigma) #MASS package


z1=1+0.5*x1+x2+e[,1]
y1=1*(z1>=0)

z2=0.8+0.3*x1+1.2*x2+e[,2]
y2=1*(z2>=0)

我试图估计的参数是潜在效用函数 z1 和 z2 中的 beta，以及方差-协方差矩阵中的非对角线元素。

谢谢！

我首先指定错误，然后在错误之后提供代码：

首先，似乎源于代码中这一行的错误：

mle = optim(theta.start,logl,x=x,y1=y1,y2=y2,hessian=T)  #Error@Here.

A）如果我在调用 OPTIM 的参数中设置 hessian = F，我会收到以下错误和回溯：

Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x),  : 
  'data' must be of a vector type, was 'NULL' 

6 array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), 
    NULL) else NULL) 
5 as.matrix.default(a) 
4 as.matrix(a) 
3 solve.default(mle$hessian) 
2 solve(mle$hessian) 
1 mle.reg(fmla, bvprobitdata)

B）如果我在调用 OPTIM 的参数中设置 hessian = T，我会收到以下错误和回溯：

Error in solve.default(mle$hessian) : 
  Lapack routine dgesv: system is exactly singular: U[1,1] = 0

3 solve.default(mle$hessian) 
2 solve(mle$hessian) 
1 mle.reg(fmla, bvprobitdata)

现在代码：

 # MLE Estimation of Bivariate Normal with correlation.

    require(Formula)
require(pbivnorm)
#Get probit data
bvprobitdata <- read.csv("/Users/...../yhbi_probitdata.csv", header = TRUE)
head(bvprobitdata,10)

#Bivariate Normal Estimation using MLE

mle.reg = function(fmla,data) {

   # Define the negative log likelihood function
    logl <- function(theta,x,y1,y2){
    y1 <- y1
    y2 <- y2
    x <- x

    #Id <- rep(1,1000)
    #x <- as.matrix(cbind(Id,x1,x2))
    beta1 <- matrix(theta[1:3],3,1)
    beta2 <- matrix(theta[4:6],3,1)
    ro <- theta[7]

    # Calculate CDFs
    temp1 <- as.matrix(cbind((x%*%beta1),(x%*%beta2))) # Create a matrix of the two cross products
    bvCDF <- pbivnorm(temp1,rho=ro) # Bivariate CDF
    xb1CDF <- pnorm(x%*%beta1)   
    Negxb1CDF <- pnorm(-(x%*%beta1))

    # Calculate Log Likelihood  - Temporarily commented out to focus debugging error in Hessian in OPTIM.
    #llik <- y1*y2*bvCDF + y1(1-y2)*log(xb1CDF-bvCDF) + (1-y1)*log(Negxb1CDF) #Calc log likelihood    
    #loglik <- sum(llik) # Sum up the log likelihoods for each observation.    
    #return(-loglik) # -ve Since OPTIM minimizes and we want to maximize loglikelihood.

    return(100)
    }


  # Prepare the data
  fml <- model.frame(fmla, data =data)
  fml
  outcome1 = rownames(attr(terms(fmla),"factors"))[1]
  outcome2 = rownames(attr(terms(fmla),"factors"))[2]
  head(data,10)
  print(outcome2)
  dfrTmp = model.frame(data)

  y1 = as.numeric(as.matrix(data[,match(outcome1,colnames(data))]))
  y2 = as.numeric(as.matrix(data[,match(outcome2,colnames(data))]))
  x = as.matrix(model.matrix(fmla, data=dfrTmp))

  # Define initial values for the parameters
  theta.start = cbind(1,1,1,1,1,1,0.5)
  # Assign names to the parameters
  names(theta.start)[1] = "b10"
  names(theta.start)[2] = "b11"
  names(theta.start)[3] = "b12"
  names(theta.start)[4] = "b20"
  names(theta.start)[5] = "b21"
  names(theta.start)[6] = "b22"
  names(theta.start)[7] = "ro"


  # Calculate the maximum likelihood
  mle = optim(theta.start,logl,x=x,y1=y1,y2=y2,hessian=T)  #Error@Here.  
  out = list(beta=mle$par,vcov=solve(mle$hessian),ll=2*mle$value)
}

print("before call")
fmla <- Formula(y1 | y2 ~x1+x2) #Create model formula
mlebvprobit = mle.reg(fmla,bvprobitdata) #Estimate coefficients for probit 
print("after call")
mlebvprobit
4

2 回答 2

0

您是否尝试过分析可能性？通常，当没有足够的参数数量数据时会发生这种情况，因此无法反转或 cholesky 变换粗麻布。

于 2014-03-26T13:57:08.677 回答
0

这有点延迟，但您想取负粗麻布的对角线的平方根，因为您的 Fisher 信息矩阵将为您提供模型的 SSE。

sqrt(solve(-mle$hessian))
于 2016-03-17T14:23:44.743 回答