[R-br] Superfície de Resposta
Emerson Cotta Bodevan
bodevan.ec em gmail.com
Terça Agosto 23 15:46:24 BRT 2016
Prezados Edson e Felipe,
obrigado pela pronta ajuda.
Att.,
*Emerson*
Em 23 de agosto de 2016 13:57, Felipe <felipe.e.barletta em gmail.com>
escreveu:
> Segue um exemplo simples com dados fictícios:
> ##### Exemplo com x1:Tempo e x2:temperatura e y:reposta
> x1 <- c(-1,-1,1,1,0,0,0,0,0,sqrt(2),-sqrt(2),0,0)
> x2 <- c(-1,1,-1,1,0,0,0,0,0,0,0,sqrt(2),-sqrt(2))
> y2 <- c(76.5,77,78,79.5,79.9,80.3,80,79.7,79.8,78.4,75.6,78.5,77)
> mod <- lm(y2 ~ x1*x2+I(x1^2)+I(x2^2))
>
> summary(mod)
>
> alr3::pureErrorAnova(mod)
>
>
> ### Para construir a superficie de resposta, inicialmente, deve ser gerada
> ### uma sequencia valores dos niveis dos fatores
> x <- seq(-sqrt(2),sqrt(2),l=30)
> y <- seq(-sqrt(2),sqrt(2),l=30)
>
> z <- function(x=x,y=y){
> y <- mod$coef[1]+mod$coef[2]*x+mod$coef[3]*y+mod$coef[4]*x^2+
> mod$coef[5]*y^2+mod$coef[6]*x*y
> }
> mod
>
> ## Comando outer prepara os dados para serem utilizados pela função persp
> z <- outer(x,y,z)
> z
> persp(x,y,z,theta=-35,phi=5,expand=0.5,xlab='Tempo',ticktype='detailed',
> ylab='Temperatura',zlab='Resposta',scale=T,col='gray80')
>
>
> image(x,y,z,col=heat.colors(155))
> contour(x,y,z,add=T)
>
> # Incluindo legenda
> fields::image.plot(x,y,z,col=rev(heat.colors(155)))
>
>
>
> ######### Ponto estacionario
> b <- matrix(c(mod$coef[2],mod$coef[3]));b
> B <- matrix(c(mod$coef[4],mod$coef[6]/2,mod$coef[6]/2,mod$coef[5]
> ),ncol=2);B
>
> x0 <- -.5*solve(B)%*%b;x0
>
> tempo <- x0[1]*5+85
> temperatura <- x0[2]*5+175
>
> image(x,y,z,col=heat.colors(115))
> contour(x,y,z,add=T)
> points(x0[1],x0[2],col='darkgreen',pch=19)
>
> # Ponto de maximo
> y0 <- mod$coef[1]+1/2%*%t(x0)%*%b
> y0
>
> --
> Atenciosamente
> Felipe E. Barletta Mendes
> Estatístico(UFPR) - Conre3 9766-A
> Mestrando em Bioestatística(UEM)
> +55 (41)-92077191
> +55 (41)-33287216
>
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160823/7779592b/attachment.html>
Mais detalhes sobre a lista de discussão R-br