Archive

Posts Tagged ‘derivada’

Gamma não central (ncg v 0.1.0)

A distribuição gama incompleta não-central pode ser vista como a forma geral da distribuição \chi^2 não-central e é expressa como uma mistura de uma Poisson com a função gamma (\Gamma) incompleta.

Alguns autores propuseram aproximações para a obtenção da função gamma incompleta não-central, sem, no entanto, que os procedimentos numéricos tivessem sido implementados. Outros trabalhos apresentam métodos para a obtenção da função de distribuição gamma (\Gamma) incompleta não-central e de sua inversa, mas os códigos de suas soluções não estão disponíveis.

Para o caso particular da distribuição \chi^2 não-central, existem soluções numéricas de alta qualidade e implementações em programas como o R (família de funções *chisq(), usando-se o argumento ncp). Entretanto, um dos problemas ainda remanescentes é o da obtenção do parâmetro de não-centralidade, dados os valores da função de distribuição, de x e de \alpha, tanto para a distribuição gama (\Gamma), quanto para a \chi^2 não-central.

Neste post apresentamos o pacote (ncg) que implementa uma adaptação ao método originalmente proposto para a função \beta incompleta não-central e um outro que combina o método desses autores com o método da inversão da função de distribuição em relação ao parâmetro de não-centralidade, utilizando para isso o método de Newton-Raphson.

Maiores detalhes acerca da teoria e eficácia do implementado no pacote estão disponíveis na dissertação que serviu de cerne ao desenvolvimento do Pacote, de autoria de Izabela Regina Cardoso de Oliveira, sob orientação do Prof Dr. Daniel Furtado Ferreira, do programa de Pós-Graduação em Estatística e Experimentação Agrícola da Universidade Federal de Lavras (texto completo aqui). Detalhes do pacote, estão também disponíveis no CRAN, (acesse-o por aqui)

Para exemplificar, suponha uma v.a. Y com distribuição \chi^2 com \nu = 5 e parâmetro de não-centralidade \delta = 1. Para calcular a probabilidade acumulada em y = 4 usa-se a função pchisq(), como segue:

## install.packages('ncg')
require(ncg)

nu <- 5
y <- 4
delta <- 1

pchisq(y, nu, 1)

Como a \chi^2 não-central é um caso particular da distribuição \Gamma não-central, por meio de uma simples transformação podemos calcular a probabilidade acumulada da \chi^2, usando a função pgammanc() do pacote ncg, fazemos portanto:

alpha <- nu / 2
x <- y / 2

pgammanc(x, alpha, delta) 

A função deltagammanc() computa o parâmetro de não-centralidade da função \Gamma não-central dados os valores de x, \alpha e p. Como a \chi^2 não-central é um caso particular da \Gamma não-central, essa função também pode ser usada para obter o parâmetro de não-centralidade de uma \chi^2 não-central.

deltagammanc(x, alpha, 0.3471897)

Espero que tanto o pacote, quanto as breves explicações aqui ditas sejam úteis. Agradeço aos colegas que estiveram envolvidos, Izabela e Prof Daniel, obrigado!

Newton-Raphson em uma dimensão

Trajetória obtida com os valores consecutivos nas iterações do método Newton-Raphson para estimação da raíz de uma função.

O método Newton-Raphson é empregado para encontrar o zero de funções. Em problemas de otimização, é empregado para encontrar o zero da função derivada f'(x) de nossa função objetivo f(x). O procedimento se baseia na aproximação em série de Taylor de primeira ordem ao redor de um valor arbitrário x_0. Considerando que desejamos obter o zero de uma função f(x), temos a expansão ao redor de x_0

f(x) \approx f(x_0) + f'(x_0)(x-x_0)

como queremos saber qual o valor de x para o qual f(x)=0, nós isolamos x na expressão acima, assim

x = x_0 + f(x_0)/f'(x_0)

e como esse valor é baseado em uma aproximação linear ao redor de x_0, agora atualizamos x_0 com o valor encontrado e repetimos esse procedimento até alcançar convergência por algum critério. Normalmente esses critérios são sobre as diferenças entre os valores obtidos em passos consecutivos.

No CMR abaixo eu apresento a obtenção raízes de uma função composta por \sin(x) e x^2. Nos escrevemos as funções, obtemos a derivada e aplicamos o método Newton-Raphson observando graficamente o histórico de iterações. No final eu coloquei um exemplo de como obter o ponto de máximo de uma função aplicando o método Newton-Raphson. Essas funções são didáticas e o método Newton-Raphson está disponível na função micEcon::maxNR() para maximização baseada em Newton-Raphson. No pacote animation a função newton.method() ilustra o funcionamento do método. Mais sobre métodos de otimização estão disponíveis na Task View Optimization do R. Em um post futuro vou ilustrar o uso do método Newton-Raphson em estimação por verossimilhança. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# newton-raphson em uma dimensão
# http://condicaoinicial.com/2011/08/sistema-nao-lineares-newton-raphson.html
# http://www.mat.ufrgs.br/~guidi/grad/MAT01169/laminas.calculo_numerico.5.pdf
# http://www.editora.ufla.br/site/_adm/upload/revista/23-3-1999_25.pdf
# http://www.dpi.inpe.br/~miguel/cursoemilioribas/MaterialRef/Curso1Estatistica_Corina.pdf

f <- quote(sin(x)-x^2/16)        # expressão da função, queremos obter a raíz
fx0 <- function(x){ eval(f) }    # função
curve(fx0, -10, 10); abline(h=0) # gráfico da função

f1 <- D(f,"x")                   # expressão da derivada
fx1 <- function(x){ eval(f1) }   # função
par(new=TRUE); curve(fx1, -10, 10, col=2, axes=FALSE) # gráfico

#-----------------------------------------------------------------------------
# método iterativo de Newton-Raphson

i <- 1       # valor inicial para o passo
dif <- 100   # valor inical para a diferença entre duas avaliações
tol <- 10^-9 # diferência mínima entre duas avaliações (tolerância)
dif <- 100   # número máximo de passos
x <- -7      # valor inicial para a raiz

while(i<100 & dif>tol){
  x[i+1] <- x[i]-fx0(x[i])/fx1(x[i])
  dif <- abs(x[i+1]-x[i])
  i <- i+1
  print(x[i])
}

#-----------------------------------------------------------------------------
# gráfico que ilustra a trajetória até a obtenção da raíz

curve(fx0, -10, 10)
for(j in 2:i){
  arrows(x[j-1], fx0(x[j-1]), x[j], fx0(x[j]), length=0.1, col=3)
  Sys.sleep(0.5)
}
abline(v=x[i], h=fx0(x[i]), col=2)

# experimente outros valores iniciais
# veja que essa função tem 2 raízes, cada valor inicial vai levar
# à apenas uma delas!

#-----------------------------------------------------------------------------
# fazendo animação gif

dir.create("seqpngs")
setwd("seqpngs")

png(file="p%02d.png", width=400, height=300)
for(j in 2:i){
  curve(fx0, -10, 10, main=paste("passo ", j-1, ", (x = ",
                        round(x[j],2), ")", sep=""))
  arrows(x[j-1], fx0(x[j-1]), x[j], fx0(x[j]), length=0.1, col=3, lwd=2)
  abline(v=x[j], h=fx0(x[j]), col="gray50")
}
abline(v=x[i], h=fx0(x[i]), col=2, lwd=2)
text(0, -3, label="CONVERGIU!", cex=2)
dev.off()

system("convert -delay 80 *.png example_1.gif")

#-----------------------------------------------------------------------------
# encontrando o ponto de máximo de uma função com newton-raphson

f <- quote(x*(2+0.5*x)^(-4))     # expressão da função, queremos obter MÁXIMO
fx0 <- function(x){ eval(f) }    # função
curve(fx0, 0, 10)

f1 <- D(f,"x")                   # expressão da derivada, obteremos o zero
fx1 <- function(x){ eval(f1) }   # função
curve(fx1, 0, 10, col=2)         # gráfico

f2 <- D(f1, "x")
fx2 <- function(x){ eval(f2) }   # função

#-----------------------------------------------------------------------------
# aplicando o newton-raphson

i <- 1       # valor inicial para o passo
dif <- 100   # valor inical para a diferença entre duas avaliações
tol <- 10^-9 # diferência mínima entre duas avaliações (tolerância)
dif <- 100   # número máximo de passos
x <- 0       # valor inicial para a raiz

while(i<100 & dif>tol){
  x[i+1] <- x[i]-fx1(x[i])/fx2(x[i])
  dif <- abs(x[i+1]-x[i])
  i <- i+1
  print(x[i])
}

#-----------------------------------------------------------------------------
# trajetória

curve(fx0, 0, 10)
for(j in 2:i){
  arrows(x[j-1], fx0(x[j-1]), x[j], fx0(x[j]), length=0.1, col=3)
  Sys.sleep(0.5)
}
abline(v=x[i], h=fx0(x[i]), col=2)

# experimente outros valores iniciais, como por exemplo 3 e 8
# para essa última função podemos obter o ponto de máximo analiticamente
#-----------------------------------------------------------------------------

Agradecimentos aos meus alunos Estevão Prado e Fillipe Pierin (acadêmicos do Curso de Estatística – UFPR) pelo desenvolvimento do núcleo código e pela permissão para publicar no ridículas.

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)
}
#-----------------------------------------------------------------------------

Gráficos de funções e suas derivadas

Nessa dica vou apresentar o CMR (código mínimo reproduzível) para fazer gráfico de funções, obter derivadas simbólicas de funções, adicionar o gráfico das derivadas das funções com escala no eixo das abcissas da direita com as expressões na legenda do gráfico (ver figura abaixo).

Gráfico da função e sua primeira derivada.

 

A curve() é a função do R que faz gráficos no plano cartesiano de função paramétrica de uma variável. Diversas opções estão disponíveis para essa função, como os tipos de linhas, cores e demais parâmetros gráficos.

A D() é a função do R que retorna a expressão das derivadas (simbólicas) de uma função. O resultado desta e um objeto de classe language. Essa função funciona de maneira recursiva como você pode observar ao obter a derivada de segunda ordem. Infelizmente, o resultado da derivação não é simplificado matematicamente.

A sobreposição de gráficos é feita nas linhas [16,25,29] e a adição de uma escala para o eixo das abcissas foi colocada no lado direito do gráfico com o comando das linhas [18,28].

A função do.call() foi usada para tornar o procedimento mais automático no sentido de que com um objeto especificando a função obteremos todos os mesmos resultados. Saber automatizar os processos é importante quando tem que e avaliar diversas funções. Imagina ficar dando copia e cola o tempo todo? Vamos deixar isso para pessoas que usam programas (em geral pagos) que não oferecem os recursos que temos com o R. Tempo é dinheiro!

Com a função uniroot() obtemos as raízes da nossa função e representamos no gráfico com uma reta vertical.

Finalmente, confeccionamos uma figura para a publicação com rótulos em todos os eixos e legenda dentro do gráfico com a expressão de cada função.

#-----------------------------------------------------------------------------
# gráfico da função, um polinômio de ordem 3 com limites em x de -10 a 10

curve(1+1/2*x+1/3*x^2+1/4*x^3, -10, 10)

#-----------------------------------------------------------------------------
# derivada de primeira e segunda order dessa expressão com relação à x

D(expression(1+1/2*x+1/3*x^2+1/4*x^3), "x")
D(D(expression(1+1/2*x+1/3*x^2+1/4*x^3), "x") ,"x")

#-----------------------------------------------------------------------------
# sobre posição do gráfico da derivada primeira adicionando um novo eixo y

curve(1+1/2*x+1/3*x^2+1/4*x^3, -10, 10)
par(new=TRUE)
curve(1/2+1/3*(2*x)+1/4*(3*x^2), -10, 10, xlab="", ylab="", col=2, axes=FALSE)
axis(4, col=2)

#-----------------------------------------------------------------------------
# procedimento mais automático, atribuindo a expressão a um objeto

expr <- expression(1+1/2*x+1/3*x^2+1/4*x^3)
do.call(curve, list(parse(text=as.character(expr)), -10, 10, ylab=expr))
par(new=TRUE)
do.call(curve, list(D(expr, "x"), -10, 10,
                    xlab="", ylab="", col=2, axes=FALSE))
axis(4, col=2)
par(new=TRUE)
do.call(curve, list(D(D(expr, "x"), "x"), -10, 10,
                    xlab="", ylab="", col=3, axes=FALSE))

#-----------------------------------------------------------------------------
# adicionando a linha da raízes do polinômio

ro <- uniroot(function(x) 1+1/2*x+1/3*x^2+1/4*x^3, c(-10, 10))
str(ro)

curve(1+1/2*x+1/3*x^2+1/4*x^3, -10, 10)
abline(v=ro$root, col=1, h=0, lty=3)

#-----------------------------------------------------------------------------
# confeccionando uma figura para publicação, para salvar descomente

#png("f002.png", w=500, h=300);
par(mar=c(5.1,4.1,2.1,4.1))
expr <- expression(1+1/2*x+1/3*x^2+1/4*x^3)
do.call(curve, list(parse(text=as.character(expr)), -10, 10,
                    ylab=expression(f(x))))
par(new=TRUE)
do.call(curve, list(D(expr, "x"), -10, 10,
                    xlab="", ylab="", lty=2, axes=FALSE))
axis(4)
mtext(side=4, line=3, text=expression(df(x)/dx))
legend("topleft", bty="n", lty=1:2,
       legend=c(paste(expression(f(x)), "=", expr),
         paste(expression(df(x)/dx), "=", capture.output(D(expr, "x")))))
#dev.off()

#-----------------------------------------------------------------------------

Veja a documentação dessas funções para um melhor aproveitamento. Deixe suas sugestões/dúvidas nos comentários. Até a próxima Ridícula!