Фитинг пробит модель INR R - PullRequest
Новая
       40

Фитинг пробит модель INR R

0 голосов
/ 05 марта 2019

Для моей дипломной работы мне нужно дополнить некоторые модели glm MLE, которых нет у R, я подходил для моделей с близкой формой, но теперь мне нужно использовать 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.Спасибо !!

1 Ответ

0 голосов
/ 05 марта 2019

Строка кода, ответственная за вашу ошибку, выглядит следующим образом: 

eta <- X %*% beta

Обратите внимание, что "% *%" является оператором умножения матриц.Воспроизводя ваш код, я заметил, что X - это матрица с 50 строками и 4 столбцами.Следовательно, чтобы умножение матриц было возможным, ваша «бета» должна иметь 4 строки.Но когда вы запускаете «betas [-p]», вы поднаборываете вектор бета-версии, удаляя последний элемент, оставляя только три элемента вместо четырех, необходимых для определения умножения матриц.Если вы удалите [-p], код будет работать.

...