optim（）を使ってプロビット回帰モデルを推定する

Question

glmを使わずにプロビット回帰モデルを手動でプログラムする必要があります。負の対数尤度を直接最小化するためにoptimを使用します。

私は以下のコードを書きましたが、それはエラーを与えて、動作しません：

cannot coerce type 'closure' to vector of type 'double'

# load data: data provided via the bottom link 
Datospregunta2a <- read.dta("problema2_1.dta") 
attach(Datospregunta2a) 

# model matrix `X` and response `Y` 
X <- cbind(1, associate_professor, full_professor, emeritus_professor, other_rank) 
Y <- volunteer 

# number of regression coefficients 
K <- ncol(X) 

# initial guess on coefficients 
vi <- lm(volunteer ~ associate_professor, full_professor, emeritus_professor, other_rank)$coefficients 

# negative log-likelihood 
probit.nll <- function (beta) { 
    exb <- exp(X%*%beta) 
    prob<- rnorm(exb) 
    logexb <- log(prob) 
    y0 <- (1-y) 
    logexb0 <- log(1-prob) 
    yt <- t(y) 
    y0t <- t(y0) 
    -sum(yt%*%logexb + y0t%*%logexb0) 
    } 

# gradient 
probit.gr <- function (beta) { 
    grad <- numeric(K) 
    exb <- exp(X%*%beta) 
    prob <- rnorm(exb) 
    for (k in 1:K) grad[k] <- sum(X[,k]*(y - prob)) 
    return(-grad) 
    } 

# direct minimization 
fit <- optim(vi, probit.nll, gr = probit.gr, method = "BFGS", hessian = TRUE) 

データ：https://drive.google.com/file/d/0B06Id6VJyeb5OTFjbHVHUE42THc/view?usp=sharing

Answer 1

大文字と小文字を区別

Yとyは異なっています。したがって、定義された関数probit.nllとprobit.grには、Yではなくyを使用する必要があります。

これらの2つの機能も私には似ていません。最も明らかな問題は、rnormの存在です。以下は正しいものです。これは次のように聞こえないlm()

から

負の対数尤度関数

# requires model matrix `X` and binary response `Y` 
probit.nll <- function (beta) { 
    # linear predictor 
    eta <- X %*% beta 
    # probability 
    p <- pnorm(eta) 
    # negative log-likelihood 
    -sum((1 - Y) * log(1 - p) + Y * log(p)) 
    } 

勾配関数

# requires model matrix `X` and binary response `Y` 
probit.gr <- function (beta) { 
    # linear predictor 
    eta <- X %*% beta 
    # probability 
    p <- pnorm(eta) 
    # chain rule 
    u <- dnorm(eta) * (Y - p)/(p * (1 - p)) 
    # gradient 
    -crossprod(X, u) 
    } 

初期パラメータ値合理的な考え。バイナリデータに線形回帰を適用する必要はありません。

しかし、純粋にlmの使用に焦点を当てると、式の右辺の共変量を分離するには+ではなく,が必要です。

---

再現性の例

彼らは互いに非常に接近している

set.seed(0) 
# model matrix 
X <- cbind(1, matrix(runif(300, -2, 1), 100)) 
# coefficients 
b <- runif(4) 
# response 
Y <- rbinom(100, 1, pnorm(X %*% b)) 

# `glm` estimate 
GLM <- glm(Y ~ X - 1, family = binomial(link = "probit")) 

# our own estimation via `optim` 
# I am using `b` as initial parameter values (being lazy) 
fit <- optim(b, probit.nll, gr = probit.gr, method = "BFGS", hessian = TRUE) 

# comparison 
unname(coef(GLM)) 
# 0.62183195 0.38971121 0.06321124 0.44199523 

fit$par 
# 0.62183540 0.38971287 0.06321318 0.44199659 

のは、おもちゃのデータセットを生成してみましょう！