Saudações. Está previsto para a disciplina de estudos dirigidos em estatística a abordagem do tópico Modelos Não Lineares. Para início de discussão os seguintes materiais devem ser lidos:
#----------------------------------------------------------- # potássio liberado acumulado para o esterco de codorna klib <- c(51.03, 57.76, 26.60, 60.65, 87.07, 64.67, 91.28, 105.22, 72.74, 81.88, 97.62, 90.14, 89.88, 113.22, 90.91, 115.39, 112.63, 87.51, 104.69, 120.58, 114.32, 130.07, 117.65, 111.69, 128.54, 126.88, 127.00, 134.17, 149.66, 118.25, 132.67, 154.48, 129.11, 151.83, 147.66, 127.30) #----------------------------------------------------------- # tempo em que foram feitas as coletas tempo <- rep(c(15, 30, 45, 60, 75, 90, 120, 150, 180, 210, 240, 270), each=3) liber <- data.frame(tempo, klib) require(lattice) #----------------------------------------------------------- # previa gráfica dos dados xyplot(klib~tempo, data=liber, type=c("p", "smooth"), col=1, xlab="Período de incubação (dias)", ylab="Potássio liberado acumulado (mg/kg de solo)") #----------------------------------------------------------- # conjunto de dados com as médias das repetições e prévia lmedio <- data.frame(tempo=unique(liber$tempo), kmedio=tapply(liber$klib, liber$tempo, mean)) xyplot(kmedio~tempo, data=lmedio, type=c("p", "smooth"), col=1, xlab="Período de incubação (dias)", ylab="Potássio liberado acumulado (mg/kg de solo)") #----------------------------------------------------------- # daqui em diante é com você... #----------------------------------------------------------- # exemplo de uso da optim() x <- 1:9 A <- 5 B <- 1 y <- A*x/(B+x)+rnorm(x,0,0.1) plot(y~x) curve(A*x/(B+x), add=TRUE) #----------------------------------------------------------- # definição da função objetivo fun.objetivo <- function(theta, y, x){ sum((y-theta[1]*x/(theta[2]+x))^2) } #----------------------------------------------------------- # escolha de valores iniciais start <- c(3,0.5) #----------------------------------------------------------- # optimização da função objetivo opt <- optim(start, fun.objetivo, y=y, x=x) opt curve(opt$par[1]*x/(opt$par[2]+x), add=TRUE, col=2) #----------------------------------------------------------- # usando outra função objetivo fun.objetivo <- function(theta, y, x){ n <- length(y) -(-n/2*log(2*pi)-n/2*log(theta[3])- sum((y-theta[1]*x/(theta[2]+x))^2/(2*theta[3]))) } #----------------------------------------------------------- # os chutes start <- c(3,0.5,0.1) #----------------------------------------------------------- # optimização opt <- optim(start, fun.objetivo, y=y, x=x) opt curve(opt$par[1]*x/(opt$par[2]+x), add=TRUE, col=3) #----------------------------------------------------------- #------------------------------------------------------------------------------------------ library(gWidgetsRGtk2) da <- data.frame(x=1:20) da$y <- 10*da$x/(3+da$x)+rnorm(da$x,0,0.2) plot(y~x, data=da) #------------------------------------------------------------------------------------------ limits <- list(A=c(0,20), B=c(0,6), n=c(0,2)) plotMM <- function(...){ plot(y~x, data=da) curve(svalue(A)*x^svalue(n)/(svalue(B)+x), add=TRUE) } w <- gwindow("Slider and spinbox example") tbl = glayout(cont=w) for(i in 1:length(limits)){ tbl[i,1] <- paste("Slide to adjuste parameter", names(limits)[i]) tbl[i,2, expand=TRUE] <- (assign(names(limits)[i], gslider(from=limits[[i]][1], to=limits[[i]][2], by=diff(limits[[i]])/20, value=mean(limits[[i]]), container=tbl, handler=plotMM))) } plotMM() #------------------------------------------------------------------------------------------