Scoring, Hessian and Newton Raphson function of the exponential distribution in R

Asked

Viewed 437 times

1

How to find the score vector and the Hessian matrix in R to apply the Newton Raphson in the code below:

rm(list=ls())
cat("\014")

#Função para simular variáveis aleatórias de uma modelo de regressão exponencial.

simula.exponencial <- function(formula, beta) {
  X <- model.matrix(formula)
  lambda <- exp(-X %*% beta)
  y <- rexp(nrow(X), lambda)
  return(data.frame(y = y, X))
}

set.seed(123)
n=10
cov <- seq(0, 5, length=n)
dados1 <- simula.exponencial(~cov, c(2,0.5))
dados1

## Função escore
escore <- function(par, formula, dados){
  mf <- model.frame(formula, dados)
  X <- model.matrix(formula, data=mf)
  esco <- ?????????
  return(drop(esco))
  }

## Hessiana

hessiano <- function(par, formula, dados){
  X <- model.matrix(formula, data=dados)
  mat <- matrix(0, nrow(X), nrow(X))
  diag(mat) <- ???????
  H <- ??????????
  return(H)
}

Thanks in advance!

r regression

  • I’m sorry, I didn’t quite understand your question. It has to do with ????? Another thing: could you provide a smaller example to demonstrate your doubt? The smaller the program, the easier it is for someone to answer the question. Even better if your program doesn’t even need the scroll bar :)

  • Okay hugomg, my question is about the parts with ??? yes, I couldn’t build these functions (Score and Hessian). I put the whole program because I believe it is interesting to help someone with this same problem! But I agree with you and I apologize for the size of the post!

  • 1

    @fsbmat, your question is very bad the way it is. This problem is not just finding a code, it has a theoretical aspect. So, my suggestions are 1) that you use Latex and define exactly the model you are sampling, to make it clear that you are not generating the data inappropriately and 2) specify theoretically what you want to find, using Latex if necessary. Note that if it is an exponential regression I could apply the log and use the maximum likelihood estimators. Downvote waiting for changes.

1 answer

0


Hi @Flaviobarros, thank you for your feedback! I managed to resolve this issue, I apologize for the delay in responding. The reply script is:

## Função escore
escore <- function(par, formula, dados){
mf <- model.frame(formula, dados)
X <- model.matrix(formula, data=mf)
esco <- crossprod((n*(t(X)*par)^(-1)-model.response(mf)), t(X))
return(drop(esco))
}

## Hessiana ("naïve")
hessiano <- function(par, formula, dados){
X <- model.matrix(formula, data=dados)
mat <- matrix(0, nrow(X), nrow(X))
diag(mat) <- -exp(X%*%par)
H <- n/par^2
return(H)
}

The same can be found in the SINAPE PDF of Prof. Paulo Justiniano and others: http://www.leg.ufpr.br/lib/exe/fetch.php/pessoais:mcie-sinape-v12.pdf

Browser other questions tagged r regression

You are not signed in. Login or sign up in order to post.