Início > inferência, regressão > Ilustração do método Newton-Raphson em regressão não linear

Ilustração do método Newton-Raphson em regressão não linear

Trajetória obtida com os valores consecutivos nas iterações do método Newton-Raphson para estimação de parâmetros em um modelo de regressão não linear.

Em regressão não linear, as estimativas dos parâmetros são obtidas por emprego de algum método numérico. Um dos métodos empregados é o de Newton-Raphson. Tal método é usado para encontrar raízes de uma função. Portanto, não é um método de maximação/mínimização de funções, mas usa o fato equivalente de que a derivada no ponto crítico de uma função é zero.

Documentação sobre o método Newton-Raphson existe em abundância. O que motivou fazer esse post foi não ter encontrado nenhum que exemplificasse aplicação para estimação de parâmetros em um modelo de regressão não linear. O estimador de mínimos quadrados para um modelo de regressão é aquele que torna mínima a função

S(\theta) = \sum_{i=1}^n (y_i-f(x_i, \theta))^2.

O método Newton-Raphson vai procurar as raízes da derivada dessa função com relação à cada parâmetro. Então, se \theta é um vetor de p parâmetros, a derivada de S(\theta) será uma função p-dimensional,

S_1(\theta) = \displaystyle\frac{\partial S(\theta)}{\partial\theta}=\begin{bmatrix}  \displaystyle\frac{\partial S(\theta)}{\partial\theta_1}\\  \vdots\\  \displaystyle\frac{\partial S(\theta)}{\partial\theta_p}  \end{bmatrix}.

O método se baseia na aproximação de primeira ordem da função S_1(\theta) de forma que um valor da função para um valor \theta^{i+1} na redondeza de um valor \theta^i é

S_1(\theta^{i+1}) \approx S_1(\theta^i)+ (\theta^{i+1}-\theta^i)\displaystyle\frac{\partial S_1(\theta)}{\partial\theta}|_{\theta=\theta^i}.

O valor de S_1(\theta^{i+1}) é zero uma vez que \theta^{i+1} é a raíz da função considerando a expansão de primeira ordem, assim, reorganizando a expressão, temos que a atualização dos valores é data por

\theta^{i+1} = \theta^i- \displaystyle\frac{S_1(\theta^i)}{S_2(\theta^i)},

sendo que

S_2(\theta) = \frac{\partial S_1(\theta)}{\partial\theta}=\begin{bmatrix}  \frac{\partial S_1(\theta)}{\partial\theta_1\partial\theta_1} & \ddots\\  \ddots& \frac{\partial S_1(\theta)}{\partial\theta_p\partial\theta_p}  \end{bmatrix}

é a matriz de derivadas que tem dimensão p\times p.

Esse procedimento numérico foi implementado no script abaixo para estimação dos parâmetros do modelo Michaelis-Menten e o modelo Exponencial reparametrizado que desenvolvi na minha Dissertação de Mestrado. No início eu monto as funções separadas e no final eu faço uma função que recebe os ingredientes mínimos e faz a estimação dos parâmetros. No final do script tem um procedimento gráfico para seleção dos chutes iniciais e acompanhar a trajetória durante as iterações. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# método Newton-Raphson
# objetivo é mínimizar a soma de quadrados dos resíduos mas iremos
# encontrar as raízes da derivada dessa função pelo método de Newton-Raphson

#-----------------------------------------------------------------------------
# parâmetros, modelo e dados simulados deles

theta <- c(A=10, B=2) # parâmetros
x <- 1:12
f <- function(theta, x){
  ## função que retorna os valore preditos para theta e x fornecidos
  ## theta: vetor de parâmetros A e B do modelo Michaelis-Menten
  ## x: vetor de valores da variável independente
  haty <- theta[1]*x/(theta[2]+x)
  return(haty)
}
set.seed(222); e <- rnorm(x,0,0.1)
y <- f(theta=theta, x=x)+e
plot(y~x)
curve(f(theta, x), add=TRUE, col=2)

#-----------------------------------------------------------------------------
# função objetivo, derivada primeira e segunda com relação à theta

S0 <- function(theta, f, y, x){
  ## função que retorna a soma de quadrados dos resíduos para
  ## theta: vetor de parâmetros do modelo f
  ## f: função do modelo de regressão não linear que depende de theta e x
  ## y e x: vetores dos valores das variáveis dependente e independente
  sqr <- sum((y-f(theta, x))^2)
  return(sqr)
}

S0(theta=theta, f=f, y=y, x=x)

# derivada primeira, dimensão 1 x p
D(expression((y-A*x/(B+x))^2), "A") # derivada de S0 com relação a A
D(expression((y-A*x/(B+x))^2), "B") # derivada de S0 com relação a B

# derivada segunda, dimensão p x p
D(D(expression((y-A*x/(B+x))^2), "A"), "A") # derivada de S0 com relação a A e A
D(D(expression((y-A*x/(B+x))^2), "B"), "B") # derivada de S0 com relação a B e B
D(D(expression((y-A*x/(B+x))^2), "A"), "B") # derivada de S0 com relação a A e B
D(D(expression((y-A*x/(B+x))^2), "B"), "A") # derivada de S0 com relação a B e A

#-----------------------------------------------------------------------------
# para ter essas funções vamos usar a deriv3()

ff <- deriv3(~(y-A*x/(B+x))^2,
             c("A","B"),
             function(x, y, A, B){ NULL })

ff(x=x, y=y, theta[1], theta[2])

#-----------------------------------------------------------------------------
# mas as matrizes são 1 x p e p x p, então falta fazer a soma em n

fff <- function(x, y, A, B){
  ## função que retorna gradiente e hessiano avaliado em theta
  ## y e x: vetores dos valores das variáveis dependente e independente
  ## A e B: valores dos parâmetros do modelo Michaelis-Menten
  m <- ff(x, y, A, B)
  g <- attr(m, "gradient"); G <- apply(g, 2, sum)
  h <- attr(m, "hessian"); h <- matrix(h, ncol=2*2)
  H <- apply(h, 2, sum); H <- matrix(H, 2, 2)
  return(list(G=G, H=H))
}

fff(x=x, y=y, theta[1], theta[2])

#-----------------------------------------------------------------------------
# agora é só fazer o Newton-Raphson

th0 <- c(9,2); tol <- 1e-10; niter <- 100;
dis <- 1; i <- 0
while(dis>tol & i<=niter){
  GH <- fff(x=x, y=y, th0[1], th0[2])
  G <- GH[[1]]
  iH <- solve(GH[[2]])
  th1 <- th0-iH%*%G
  dis <- sum((th0-th1)^2)
  i <- i+1
  print(c(i, dis, th1))
  th0 <- th1
}

n0 <- nls(y~A*x/(B+x), start=c(A=9, B=2))
coef(n0)

iH/vcov(n0) # hessiano é proporcional à matriz de covariância dos parâmetros

#-----------------------------------------------------------------------------
# vamos incluir todos os elementos em uma única função!

newton <- function(formula, data, start, tol=1e-10, niter=100, trace=TRUE){
  ## formula: expressão da resposta como função das covariáveis e parâmetros
  ## data: data.frame que contém dos dados
  ## start: vetor ou lista nomeados de valores iniciais para os parâmetros
  ## tol: norma mínima entre valores consecutivos dos parâmetros
  ## niter: número máximo de iterações do algortimo
  ## trace: se TRUE imprime o histórico de iterações
  form <- paste("~(",paste(formula[2:3], collapse="-"),")^2", sep="")
  vars <- all.vars(formula)
  nvars <- length(vars)
  pars <- 1:(nvars+1)
  l <- as.list(pars)
  names(l)[1:nvars] <- vars
  f <- as.function(l)
  body(f) <- NULL
  ff <- deriv3(as.formula(form), names(start), f)
  th0 <- unlist(start)
  i <- 0; dis <- 100; iter <<- matrix(th0, ncol=length(start))
  while(dis>tol & i<=niter){
    m <- do.call(ff, c(as.list(data), as.list(th0)))
    g <- attr(m, "gradient"); G <- apply(g, 2, sum)
    h <- attr(m, "hessian"); h <- matrix(h, ncol=length(start)^2)
    H <- apply(h, 2, sum); H <- matrix(H, nrow=length(start))
    iH <- solve(H)
    th1 <- th0-iH%*%G
    dis <- sum((th0-th1)^2)
    i <- i+1
    iter <<- rbind(iter, t(th1))
    if(trace==TRUE) print(c(i, dis, th1))
    th0 <- th1
  }
  th0 <- as.vector(th0)
  names(th0) <- names(start)
  ## retorna o último valor iterado e a uma matriz com o histórico
  return(list(theta=th0, iter=iter))
}

#-----------------------------------------------------------------------------
# vamos estudar vetor de chutes presentes nessa lista

start.l <- list(c(A=10,B=2.5), c(A=10,B=0.1), c(A=14,B=2.5),
                c(A=15.13, B=1.7), c(A=18.4,B=0.13))

data <- data.frame(x=x, y=y)
formula <- y~A*x/(B+x)
start <- start.l[[4]]
nr <- newton(formula, data, start)
nr$theta

sqe <- function(A, B, y, x){ hy <- (A*x)/(B+x); sum((y-hy)^2) }
SQE <- Vectorize(sqe, c("A", "B"))

A.grid <- seq(-3,15,l=100)
B.grid <- seq(-3,15,l=100)
sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
#png("f023.png", w=500, h=500)
contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
        xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
for(i in seq(nrow(iter)-1)){
  arrows(iter[i,1], iter[i,2],
         iter[i+1,1], iter[i+1,2],
         col=2, length=0.1)
  Sys.sleep(0.1)
}
#dev.off()

#-----------------------------------------------------------------------------
# com o mouse selecione regiões da superfície para serem valores iniciais

loc <- locator()
loc <- as.data.frame(loc)
names(loc) <- c("A","B")
start.l <- split(loc, f=1:nrow(loc))

#-----------------------------------------------------------------------------
# ver como fica com o modelo Exponencial

n0 <- nls(y~A*(1-exp(-log(2)*x/B)), start=c(A=9, B=2))
coef(n0)

data <- data.frame(x=x, y=y)
formula <- y~A*(1-exp(-log(2)*x/B))
start <- list(A=5, B=0.5)
nr <- newton(formula, data, start)
nr$theta

sqe <- function(A, B, y, x){ hy <- A*(1-exp(-log(2)*x/B)); sum((y-hy)^2) }
SQE <- Vectorize(sqe, c("A", "B"))

A.grid <- seq(-10,20,l=100)
B.grid <- seq(-10,20,l=100)
sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
        xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
for(i in seq(nrow(iter)-1)){
  arrows(iter[i,1], iter[i,2],
         iter[i+1,1], iter[i+1,2],
         col=2, length=0.1)
  Sys.sleep(0.1)
}
#-----------------------------------------------------------------------------

About these ads
  1. Nenhum comentário ainda.
  1. No trackbacks yet.

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s

Seguir

Obtenha todo post novo entregue na sua caixa de entrada.

Junte-se a 52 outros seguidores

%d blogueiros gostam disto: