私はあまり密度の濃いコードに反映されている可能性がある R 初心者です。ご了承ください。max を使用して二変量正規分布の係数を推定しようとしています。尤度推定。OPTIM 関数を呼び出すと、ヘッシアンに関連するエラーが発生します。かなりデバッグを試みましたが、エラーを取り除くことができないようです。これに取り組む方法についての洞察をいただければ幸いです。
私が使用しているデータは{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 のベータと、分散共分散行列の非対角要素です。
ありがとう！
最初にエラーを指定してから、エラーの後にコードを提供します。
まず、コード内のこの行で発生したと思われるエラー:
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