#============================================================================= # II CONEST - Curitiba/PR - de 7 à 9 de Novembro de 2013 # # CURSO EM MODELOS DE REGRESSÃO NÃO LINEAR # www.leg.ufpr.br/conest2013 # # Prof. DSc. Walmes Marques Zeviani # LEG - DEST - UFPR # # script04 - Uso de informação de gradiente (09/11/2013) # * Função com atributo gradiente; # * Uso da função deriv3(); # * Gradientes numéricos para modelos segmentados; # * Gradiente analítico vs gradiente numérico; #============================================================================= #----------------------------------------------------------------------------- # Definições da sessão. require(lattice) require(latticeExtra) require(segmented) require(rootSolve) require(rpanel) #----------------------------------------------------------------------------- # Usar a informação de gradientes, em geral, acelera o procedimento de # estimação. Para usar informação de gradiente função precisa ter um # atributo gradiente. data(wtloss, package="MASS") p0 <- xyplot(Weight~Days, data=wtloss, col=1, xlab="Dias em dieta", ylab="Peso (kg)") print(p0) p0+layer(with(list(f0=81, f1=102, K=141), panel.curve(f0+f1*2^(-x/K), add=TRUE, col=2))) #----------------------------------------------------------------------------- # Modelo (é uma parametrização do monomolecular para tempo de meia vida) # f0+f1*2^(-x/K) # * f0: assíntota, peso final após estabilização # * f1: total de peso à perder, f0+f1 é intercepto ou peso inicial # * K: tempo necessário para se perder metade do passível de perder, # 0.5*(f1-f0) # Para passar o gradiente tem-se que conhecer as derivadas do modelo com # relação à cada parâmetro: model <- expression(f0+f1*2^(-x/K)) D(model, "f0") D(model, "f1") D(model, "K") pars <- c("f0","f1","K") sapply(pars, D, expr=model) # todas de uma vez #----------------------------------------------------------------------------- # Agora precisa-se fazer uma função que avalie essas expressões em x e # pars e passe isso como atributo "gradient". modelfun <- function(x, f0, f1, K){ value <- f0+f1*2^(-x/K) attr(value, "gradient") <- cbind(f0=1, f1=2^(-x/K), K=f1*(2^(-x/K)*(log(2)*(x/K^2)))) return(value) } # avalia a função modelfun(x=wtloss$Days[1:5], f0=81, f1=102, K=141) #----------------------------------------------------------------------------- # Agora é só usar. # com informação de gradiente n0 <- nls(Weight~modelfun(Days, f0, f1, K), data=wtloss, start=list(f0=60, f1=200, K=41), trace=TRUE) summary(n0) # sem informação de gradiente n1 <- nls(Weight~f0+f1*2^(-Days/K), data=wtloss, start=list(f0=60, f1=200, K=41), trace=TRUE) summary(n1) # NOTA: para o presente exemplo não houve diferença em número de passos. #----------------------------------------------------------------------------- # A forma mais rápida de passar informação de gradiente é usar a função # deriv3(). Tem esse nome porque faz as derivadas do modelo com relação # aos parâmetros e retorna a avaliação de 3 funções avaliadas em x e # theta: 1) modelo (f(x)), 2) o gradiente (df(x)/dtheta) e 3) o hessiano # (d²f(x)/dtheta²). As derivadas são analíticas. Sendo assim, a deriv3() # não sabe derivar modelos com segmentados pois não sabe interpretar as # inequações do tipo x<=x0, por exemplo. modelfun <- deriv3(~f0+f1*2^(-x/K), # modelo c("f0","f1","K"), # parâmetros do modelo function(x, f0, f1, K){ NULL }) modelfun # veja que retorna 2 atributos: gradient e hessian modelfun(x=wtloss$Days[1:5], f0=81, f1=102, K=141) #----------------------------------------------------------------------------- # Alguns modelos não lineares possuem operadores/funções que não estão # na tabela de derivadas considerada pela deriv3(). Um exemplo é o # modelo linear-platô que possui os operadores < e >= na função. Dessa # forma não tem como usar informação de derivada analítica por meio da # deriv3(). Pode-se fazer à mão ou pode-se empregar derivadas # numéricas. A função rootSolve::gradiente() pode ser usada. # mensagem de erro! mlp <- deriv3(~f0+tx*x*(x=xde), c("f0","tx","xde"), function(x, f0, tx, xde){ NULL }) f <- function(theta, xi){ with(as.list(theta), f0+tx*xi*(xi=xde)) } f(c(f0=10, tx=1, xde=3), xi=1:4) # é um gradiente numérico (rootSolve) gradient(f, x=c(f0=3, tx=0.9, xde=3.5), xi=seq(0,5,by=0.5)) #----------------------------------------------------------------------------- # Criando uma função cujo gradiente é obtido numericamente. Existe # custo em obter gradiente numérico, melhor seria analítico. modelfun.n <- function(x, f0, tx, xde){ f <- function(theta, xi){ with(as.list(theta), f0+tx*xi*(xi=xde)) } value <- f(theta=c(f0=f0, tx=tx, xde=xde), xi=x) attr(value, "gradient") <- gradient(f, x=c(f0=f0, tx=tx, xde=xde), xi=x) value } modelfun.n(x=seq(0,5,by=0.5), f0=3, tx=0.9, xde=3.5) # feito! #----------------------------------------------------------------------------- # Criando uma função cujo gradiente é obtido analiticamente. modelfun.a <- function(x, f0, tx, xde){ f <- function(theta, xi){ with(as.list(theta), f0+tx*xi*(xi=xde)) } xo <- x