在R中拟合概率模型答案

【问题标题】：Fitting probit model inr R在R中拟合概率模型
【发布时间】：2019-07-26 09:09:05
【问题描述】：

对于我的论文，我必须用 R 没有的 MLE 拟合一些 glm 模型，我对形式接近的模型没问题，但现在我必须使用 de Gausian CDF，所以我决定拟合一个简单的概率模型。这是代码：

Data:
set.seed(123)
x <-matrix( rnorm(50,2,4),50,1)
m <- matrix(runif(50,2,4),50,1)
t <- matrix(rpois(50,0.5),50,1)
z <- (1+exp(-((x-mean(x)/sd(x)))))^-1 + runif(50)
y <- ifelse(z < 1.186228, 0, 1)

data1 <- as.data.frame(cbind(y,x,m,t))

myprobit <- function (formula, data) 
{
  mf <- model.frame(formula, data)
  y <- model.response(mf, "numeric")
  X <- model.matrix(formula, data = data)
  if (any(is.na(cbind(y, X)))) 
    stop("Some data are missing.")
  loglik <- function(betas, X, y, sigma) {     #loglikelihood
    p <- length(betas)
    beta <- betas[-p]
    eta <- X %*% beta
    sigma <- 1    #because of identification, sigma must be equal to 1
    G <- pnorm(y, mean = eta,sd=sigma)
    sum( y*log(G) + (1-y)*log(1-G))
  }
  ls.reg <- lm(y ~ X - 1)#starting values using ols, indicating that this model already has a constant
  start <- coef(ls.reg)


  fit <- optim(start, loglik, X = X, y = y, control = list(fnscale = -1), method = "BFGS", hessian = TRUE) #optimizar
  if (fit$convergence > 0) {
    print(fit)
    stop("optim failed to converge!") #verify convergence
  }

  return(fit)
}

myprobit(y ~ x + m + t,data = data1)

我得到：Error in X %*% beta : non-conformable arguments，如果我将start <- coef(ls.reg) 更改为start <- c(coef(ls.reg), 1)，我得到错误的提示：

probit <- glm(y ~ x + m + t,data = data1 , family = binomial(link = "probit"))

我做错了什么？是否可以使用 pnorm 正确拟合此模型，如果没有，我应该使用什么算法来近似去高斯 CDF。谢谢！！

【问题讨论】：

通常 R 人希望看到带有样本数据的完整（即：可执行）示例。见stackoverflow.com/questions/5963269/…
删除了我的“答案”，因为它显然没有响应您的代码审查类型问题。但仍然认为尚不清楚为什么 ypu 会这样做。使用能力为glm 建立新家庭有哪些障碍？
我的研究是关于罕见事件（特别是与信用相关的事物），所以我开发了一个方程系统，当你插入分布时，它们会返回一些所需的概率（国王校正是一种特殊情况，因为例如），但模型必须以某种方式拟合，对我来说，自然选择是 MLE，所以我派生了一些 MLE，现在我正在着手实现它们。但首先我必须看看我所做的是否适用于简单模型，特别是 CDF 没有紧密形式的模型。
顺便说一下，我是一个新的 R 用户（我从 3 周前开始），我将时间分配在学习 R 和拟合模型之间

标签： r statistics glm mle

【解决方案1】：

导致您的错误的代码行如下：

eta <- X %*% beta

请注意，“%*%”是矩阵乘法运算符。通过复制您的代码，我注意到 X 是一个具有 50 行和 4 列的矩阵。因此，要使矩阵乘法成为可能，您的“beta”需要有 4 行。但是当您运行“betas[-p]”时，您通过删除它的最后一个元素来对 betas 向量进行子集化，只留下三个元素而不是定义矩阵乘法所需的四个元素。如果您删除 [-p] 代码将起作用。

【讨论】：