Arquivo

Posts Tagged ‘wireframe’

Gráfico tridimensional com contornos de nível

Peso de 100 grãos em função do nível de saturação de água no solo e conteúdo de potássio fornecido pela adubação. Linhas no plano inferior representam a projeção das linhas de mesmo valor para peso de 100 grãos.

Nesse post vou apresentar uma função para adicionar contornos de nível em gráficos tridimensionais gerados pela lattice::wireframe(). Para ilustrar o uso da função eu fiz um ajuste de modelo de regressão múltipla aos dados do experimento com soja. Esses dados eu já usei em várias matérias aqui no Ridículas. Uma vez ajustado o modelo, os valores preditos são então representados pela wireframe() usando o panel.3d.contour que foi uma função que aprimorei das leituras que fiz das mensagens da r-help.

Um gráfico tridimensional como este impressiona quem o vê, porém esteja ciente das más impressões que um gráfico tridimensional pode dar pois ele não é de fato tridimensional e sim uma projeção tridimensional em um plano, o que acarreta algumas deformações. Todavia, achei por bem documentar e expor a função mas isso não significa que esta é a melhor representação. A função permite colocar os contornos no plano abaixo, acima e sobre a superfície (bottom, top, on). Tenha cuidado ao escolher o ângulo de observação do gráfico (screen=) e controle a amplitude do eixo z (zlim=) para que a superfície não esconda a projeção abaixo. Confira os resultados no CMR abaixo. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# leitura dos dados

soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/cnpaf/soja.txt",
                   header=TRUE, sep="\t", dec=",")
str(soja)

#-----------------------------------------------------------------------------
# ver o peso de 100 grãos

require(lattice)

xyplot(pesograo~potassio, groups=agua, data=soja, type=c("p","a"))
xyplot(pesograo~potassio|agua, data=soja, type=c("p","a"))

#-----------------------------------------------------------------------------
# ajuste de um modelo polinomial nos dois fatores

m0 <- lm(pesograo~bloco+poly(agua,2)*poly(potassio,2), data=soja)
par(mfrow=c(2,2)); plot(m0); layout(1)

anova(m0)
summary(m0)

m1 <- lm(pesograo~bloco+(agua+I(agua^2))*potassio+I(potassio^2), data=soja)
par(mfrow=c(2,2)); plot(m1); layout(1)

anova(m0, m1)
summary(m1)

#-----------------------------------------------------------------------------
# fazendo a predição

pred <- expand.grid(bloco="I", agua=seq(37.5,62.5,l=30),
                    potassio=seq(0,180,l=30))
pred$y <- predict(m1, newdata=pred)

#-----------------------------------------------------------------------------
# vendo a superfície ajustada

source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functions.R")
args(panel.3d.contour)
body(panel.3d.contour)

require(RColorBrewer)

colr <- brewer.pal(11, "Spectral")
colr <- colorRampPalette(colr, space="rgb")

zlab <- "Peso de 100 grãos"
xlab <- "Potássio no solo"
ylab <- "Nível de saturação de água"

#png(file="f034.png", width=500, height=500)
wireframe(y~potassio*agua, data=pred,
          scales=list(arrows=FALSE), zlab=list(zlab, rot=90),
          xlab=list(xlab, rot=24), ylab=list(ylab, rot=-37),
          zlim=c(5,16), col="gray50", col.contour=1,
          panel.3d.wireframe="panel.3d.contour", type=c("on","bottom"),
          col.regions=colr(100),  drape=TRUE,
          screen=list(z=40, x=-70))
#dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:,

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.

Como adicionar legendas em gráficos da lattice

Gráfico com bandas de confiança, legenda e assinatura de autor e versão.

No post Como fazer legendas em gráficos eu mostrei procedimentos para gráficos do pacote graphics. Nesse post vou apresentar procedimento para fazer legendas em gráficos da lattice, desenvolvido por Deepayan Sarkar. Iremos empregar as opções auto.key=, key=, colorkey=, e as funções draw.key() e draw.colorkey(). O pacote grid, desenvolvido por Paul Murrell será necessário para posicionarmos as legendas na página gráfica. Nos exemplos temos legendas para pontos, linhas, pontos e linhas, barras, escala de cores para gráficos de três dimensões, gráfico com duas legendas, posicionamento da legenda com o mouse e gráfico com legenda para representar o ajuste de um modelo de regressão com valores observados, preditos e intervalo de confiança.

#-----------------------------------------------------------------------------
# pacotes

require(lattice) # para fazer os gráficos
require(grid)    # para posicionar as legendas

#-----------------------------------------------------------------------------
# dados de indice agronômico em função da cultivar de milho e dose de N

da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                 header=TRUE, sep="\t")
str(da)

#-----------------------------------------------------------------------------
# legenda para gráfico de dispersão

# gráfico com legenda padrão
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE, auto.key=TRUE)

# gráfico com legenda em colunas e alteração dos padrões dos pontos
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
       auto.key=list(columns=3),
       par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))

# muda a posição da legenda
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
       auto.key=list(space="right", columns=1),
       par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))

# coloca a legenda dentro do gráfico
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
       auto.key=list(corner=c(0.95,0.05), columns=1),
       par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))

#-----------------------------------------------------------------------------
# legenda para gráfico de linhas

db <- with(da, aggregate(indice, list(cultivar=cultivar, dose=dose), mean))
col <- trellis.par.get("superpose.line")$col[1:nlevels(db$cultivar)]

# gráfico com legendas de linhas, nomes *antes* das linhas
xyplot(x~dose, groups=cultivar, data=db, type="l",
       key=list(corner=c(0.95,0.05), columns=1,
         text=list(levels(db$cultivar)),
         lines=list(col=col)))

# gráfico com legendas de linhas, nomes *depois* das linhas
xyplot(x~dose, groups=cultivar, data=db, type="l",
       key=list(corner=c(0.95,0.05), columns=1,
         lines=list(col=col),
         text=list(levels(db$cultivar))))

#-----------------------------------------------------------------------------
# legenda para gráfico de linhas e pontos

xyplot(x~dose, groups=cultivar, data=db, type="o", pch=19,
       key=list(corner=c(0.95,0.05), columns=1,
         type="o", divide=1,
         lines=list(col=col, pch=19),
         text=list(levels(db$cultivar))))

#-----------------------------------------------------------------------------
# legenda para pontos e reta de regressão

# ajuste de modelo e criação de objeto com valores observados e preditos
m0 <- lm(indice~cultivar*(dose+I(dose^2)), data=da)
dc <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,l=100))
dc$indice <- predict(m0, newdata=dc)
dc$ind <- "predito"
da$ind <- "obsevado"
dc <- rbind(da[,-3], dc)
str(dc)

# legenda fora do gráfico
xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("p","l"), col=1, layout=c(3,1),
       key=list(space="top", columns=2, type="o", divide=1,
         lines=list(pch=c(1,NA), lty=c(0,1)),
         text=list(c("valores observados","valores preditos"))))

# legenda dentro do gráfico
xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("p","l"), col=1, layout=c(3,1),
       key=list(corner=c(0.95,0.05), columns=1, type="o", divide=1,
         lines=list(pch=c(1,NA), lty=c(0,1)),
         text=list(c("valores observados","valores preditos"))))

#-----------------------------------------------------------------------------
# dados de mortalidade de aves em função do tipo de resfriamento do aviário

ma <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/mortes.txt",
                 header=TRUE, sep="\t")
str(ma)

#-----------------------------------------------------------------------------
# legenda para gráficos de barras

mb <- with(ma, aggregate(mortes, list(galpao=galpao, sistema=asper), sum))
str(mb)

col <- c("forestgreen", "palegreen")

# legenda fora da região gráfica
barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(space="top", columns=2,
           rectangles=list(col=col),
           text=list(c("com aspersão","sem aspersão"))))

# legenda fora da região gráfica com título
barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(title="Sistema de resfriamento", cex.title=1.1,
           space="top", columns=2,
           rectangles=list(col=col),
           text=list(c("com aspersão","sem aspersão"))))

# legenda dentro da região gráfica com título
barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(title="Sistema de resfriamento", cex.title=1.1,
           corner=c(0.05,0.95), columns=1, padding.text=2,
           rectangles=list(col=col, size=3, height=0.8, border=TRUE),
           text=list(c("com aspersão","sem aspersão"))))

#-----------------------------------------------------------------------------
# escolhendo a posição da legenda com o mouse

barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col)
trellis.focus(name="page")
loc <- grid.locator(unit="npc") # clicar no gráfico
trellis.unfocus()

barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(title="Sistema de\nresfriamento", cex.title=1.1,
           x=loc$x, y=loc$y, corner=c(0.5,0.5), columns=1, padding.text=2,
           rectangles=list(col=col, size=3, height=0.8, border=TRUE),
           text=list(c("com\naspersão","sem\naspersão"))))

#-----------------------------------------------------------------------------
# nota de rodapé no gráfico com assinatura, data e versão

# ajuste de modelo e criação de objeto com valores observados, preditos e IC
m0 <- lm(indice~cultivar*(dose+I(dose^2)), data=da)
dc <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,l=100))
indice <- predict(m0, newdata=dc, interval="confidence")
dc <- cbind(dc, stack(data.frame(indice)))
names(dc)[3] <- "indice"
dc <- rbind(da[,-3], dc)

#png("f017.png", w=500, h=350)
xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("l","l","p","l"),
       col=1, lty=c(1,2,0,2), layout=c(3,1),
       strip=strip.custom(bg="gray90"),
       xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
       key=list(space="top", columns=3, type="o", divide=1,
         lines=list(pch=c(1,NA,NA), lty=c(0,1,2)),
         text=list(c("valores observados","valores preditos","IC 95%"))))
pushViewport(viewport())
grid.text(label=paste("Walmes Zeviani --",
            format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
            R.version.string), rot=90,
          x=unit(1, "npc")-unit(2, "mm"),
          y=unit(0, "npc")+unit(2, "mm"),
          just=c("left", "bottom"),
          gp=gpar(cex=0.8, col="gray50"))
popViewport()
#dev.off()

#-----------------------------------------------------------------------------
# legenda para gráficos de 3 dimensões

x <- seq(-10,10,l=50)
dd <- expand.grid(x=x, y=x)
dd$z <- with(dd, 500+x+2*y-0.5*x*y-x^2-y^2)

# gráfico sem cores e sem legenda
wireframe(z~x+y, data=dd)

# gráfico com cores e legenda
wireframe(z~x+y, data=dd, drape=TRUE)

# remove a legenda
wireframe(z~x+y, data=dd, drape=TRUE, colorkey=FALSE)

# coloca a legenda embaixo
wireframe(z~x+y, data=dd, drape=TRUE, colorkey=list(space="bottom"))

# redimensionando a legenda
wireframe(z~x+y, data=dd, drape=TRUE,
          colorkey=list(width=1, height=0.7))

# gráfico de níveis de cor
levelplot(z~x+y, data=dd)

# com contornos e redimensiomento da legenda
levelplot(z~x+y, data=dd, contour=TRUE,
          colorkey=list(width=1, height=0.7))

#-----------------------------------------------------------------------------
# usando a draw.key(), usar fora da função gráfica ou dentro da panel()

xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("l","l","p","l"),
       col=1, lty=c(1,2,0,2), layout=c(2,2),
       strip=strip.custom(bg="gray90"),
       xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico")
draw.key(list(space="top", columns=1, type="o", divide=1,
         lines=list(pch=c(1,NA,NA), lty=c(0,1,2)),
         text=list(c("valores observados","valores preditos","IC 95%"))),
         draw=TRUE,
         vp=viewport(x=unit(0.75, "npc"), y=unit(0.65, "npc")))

#-----------------------------------------------------------------------------
# usando a draw.colorkey(), usar fora da função gráfica ou dentro da panel()

wireframe(z~x+y, data=dd, drape=TRUE, colorkey=FALSE,
          panel=function(...){
            panel.wireframe(...)
            draw.colorkey(list(space="bottom", at=seq(250,500,25),
                               height=0.8, width=1), draw=TRUE,
                          vp=viewport(x=unit(0.5,"npc"), y=unit(0.07,"npc")))
          })

#-----------------------------------------------------------------------------
# usando duas legendas

de <- expand.grid(w=gl(3,5,la="w"), z=gl(4,3,la="z"))
de$x <- with(de, rnorm(w, as.numeric(w), 0.5))
de$y <- with(de, rnorm(z, as.numeric(z), 0.5))
pch <- c(1,2,5,7); col=1:3

xy <- xyplot(y~x, data=de, groups=z:w,
             pch=rep(pch, e=3), col=rep(col, 4))

print(xy, position=c(0,0,0.8,1), more=FALSE)
draw.key(list(text=list(levels(de$z)), points=list(pch=pch, col=1),
         space="top", columns=1, title="Níveis de Z", cex.title=1.1),
         draw=TRUE, vp=viewport(x=unit(0.85, "npc"), y=unit(0.6, "npc")))
draw.key(list(text=list(levels(de$w)), points=list(pch=15, col=col),
         space="top", columns=1, title="Níveis de W", cex.title=1.1),
         draw=TRUE, vp=viewport(x=unit(0.85, "npc"), y=unit(0.4, "npc")))

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

RStudio e a função manipulate()

Manipulando ângulo de observação e gradiente de cores de um gráfico tridimensional com as ferramentas do RStudio.

O RStudio é um dos vários disponíveis editores para script R. Uma das suas peculiaridades é o pacote manipulate. Esse pacote é parte integrante do RStudio e fica disponível apenas com a sua instalação.

A função manipulate::manipulate() permite que você use os gráficos do R de forma interativa. Ou seja, através de deslizadores (slider), seletores (picker) e opções (checkbox) você é capaz de mudar aspectos gráficos. A documentação da função traz exemplos interessantes. Nessa ridícula eu implementei dois exemplos. O primeiro é para escolher o ângulo de observação de uma superfície de resposta (inclui escolha do gradiente de cores). O segundo permite que você entenda como funciona a estimação kernel de densidade.

Com a manipulate() fica muito mais simples construir funções para obter valores iniciais em modelos de regressão não linear. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados para um gráfico tridimensional

da <- expand.grid(x=seq(0,10,l=30), z=seq(0,10,l=30))
da$y <- with(da, x+z+0.2*x*z) # gera dados

#-----------------------------------------------------------------------------
# escolher o ângulo de observação e o esquema de cores
 
manipulate({
             ## faz o gráfico tridimensional
             colr <- colorRampPalette(c(c1, c2, c3), space="rgb")
             arrows <- arr
             wireframe(y~x+z, data=da, scales=list(arrows=arrows),
                       col="gray30", col.contour="gray30",
                       col.regions=colr(100),  drape=TRUE,
                       screen=list(z=z.angle, x=x.angle)
                       )
           },
           ## controla o valor dos angulos e das cores
           z.angle=slider(0, 360, step=10, initial=40),
           x.angle=slider(-180, 0, step=5, initial=-60),
           arr=checkbox(FALSE, "show.arrows"),
           c1=picker("red","yellow","orange","green","blue","pink","violet"),
           c2=picker("red","yellow","orange","green","blue","pink","violet"),
           c3=picker("red","yellow","orange","green","blue","pink","violet")
           )

#-----------------------------------------------------------------------------
# gráfico de densidade controlando o bandwidth e tipo de função kernel

x <- rgamma(300, 3, 7)

manipulate({
             plot(density(x, bw=bw, kernel=kernel))
             if(show.rug==TRUE) rug(x)
           },
           kernel=picker("gaussian", "epanechnikov", "rectangular",
             "triangular", "biweight","cosine",
             "optcosine"),
           bw=slider(0.01, 0.15, step=0.003, initial=0.05),
           show.rug=checkbox(TRUE, "show rug")
           )

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