Archive

Posts Tagged ‘density’

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!

Gráfico de densidade com destaque dos quantis

Gráficos de densidade com regiões destacando os quartis.

Gráficos de densidade são utilizados para se conhecer a forma da distribuição dos dados. A função lattice::densityplot() faz esse tipo de gráfico. Para aprimorar a visualização dos quartis dos dados, eu implementei uma função (my.densityplot()) para que destaque com regiões coloridas as porções de dados delimitadas pelos quartis. O código abaixo gera a figura no início desse post. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# densityplot() com regiões pintadas delimitadas pelos quantis

require(lattice)

my.densityplot <- function(x, cols="gray90", probs=c(1,3)/4, ...){
  panel.densityplot(x, ...)
  dx <- density(x)
  breaks <- sort(c(range(dx$x), quantile(x, probs=probs)))
  fx <- approxfun(dx$x, dx$y)
  do.polygon <- function(x){
    y <- fx(x)
    return(list(x=c(min(x), x, max(x)), y=c(0, y, 0)))
  }
  seqs <- lapply(1:(length(breaks)-1),
                 function(i){
                   x <- seq(breaks[i], breaks[i+1], l=20)
                   do.polygon(x)
                 })
  for(i in 1:length(seqs)){
    seqs[[i]]$col <- cols[i]
  }
  lapply(seqs, function(i) do.call(panel.polygon, i))
  panel.mathdensity(dmath=dnorm, col="black",
                    args=list(mean=mean(x),sd=sd(x)), lty=3)
}

require(RColorBrewer)
display.brewer.all()

cols <- brewer.pal(3, "BuPu") # Greens, Blues, BuPu, Reds
densityplot(~height|voice.part, data=singer,
            cols=cols, plot.points="rug", panel=my.densityplot)

#png("f036.png", width=500, height=320)
cols <- brewer.pal(4, "Greens") # Greens, Blues, BuPu, Reds
densityplot(~height|voice.part, data=singer, layout=c(4,2),
            strip=strip.custom(bg="gray90"),
            xlab="Altura", ylab="Densidade", col=1,
            cols=cols, probs=1:3/4, plot.points="rug", panel=my.densityplot)
#dev.off()

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

Procedimento gráfico interativo para densidade com rpanel

Printscreen da janela gráfica e janela de controle do painel para o gráfico de densidade kernel usando funções disponíveis no pacote rpanel.

Alguns posts atrás eu apresentei procedimentos gráficos para obtenção de valores iniciais em regressão não linear. Dessa vez eu mudei o foco. Vou apresentar os procedimentos gŕaficos interativos disponíveis no pacote rpanel ilustrando o procedimento de estimação kernel de densidade.

O pacote rpanel possui funções para alterar componentes do gráfico. Os mais úteis são: deslizadores, caixas de seleção e caixas de texto. Exemplos de uso dessas funções você encontra na documentação do pacote e na página do pacote mantida pelo autor do mesmo, Adrian Bowman. Para ver o rpanel funcionando, assista o vídeo abaixo.

No CMR abaixo eu preparei um painel para controlar os argumentos kernel= e width= da função density(). Além disso, outro painel controla a posição do centro da banda. Ao alterar a opção de kernel, muda-se a função usada na estimação. Ao alterar o valor da largura de banda, muda-se a zona de influência ao redor do centro. Por último, movendo o centro sobre as observações, você verifica como a função se comporta em cada região.

No final, tem uma pequena porção de código que ilustra como obter probabilidades a partir de uma função de densidade kernel. O procedimento consiste em aproximar a função kernel e depois fazer uma integração numérica.

Para usar o rpanel, basta instalar da forma usual o pacote e suas dependências. Verifique a documentação do pacote para mais opções. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados de redimento de grãos de 32 híbridos de milho

rend <- c(6388, 6166, 6047, 5889, 5823, 5513, 5202, 5172, 5166, 4975, 4778,
          4680, 4660, 5403, 5117, 5063, 4993, 4980, 4770, 4685, 4614, 4552,
          3973, 4550, 5056, 4500, 4760, 5110, 4960, 4769, 4849, 5230)

#-----------------------------------------------------------------------------
# tipos de kernel, tirado da documentação da função density()

kernels <- eval(formals(density.default)$kernel) # extrai os tipos de kernel
kernels
plot(density(0, from=-1.2, to=1.2, width=1, kernel="gaussian"), type="l",
     ylim=c(0, 2), xlab="rendimento de grãos", main="kernels")
for(i in 2:length(kernels))
  lines(density(0, width=1, kernel=kernels[i]), col=i)
legend(0.5, 2.0, legend=kernels, col=seq(kernels), lty=1, bty="n")

#-----------------------------------------------------------------------------
# painel para gráfico de densidade: largura de banda e tipo de função kernel

require(rpanel)

density.panel <- function(panel){
  ## argumento panel é uma lista, cada slot é um argumento
  ## x: vetor métrico de observações
  ## width: escalar métrico da largura de banda
  ## kernel: scalar string que informa o tipo de função kernel
  ## c0: escalar métrico que é a coodenada do centro da banda
  ## n: tamanho da amostra
  ##------------------------------------------------------------------
  ## gráfico de densidade e rug
  aux <- density(panel$x, width=panel$width, kernel=panel$kernel)
  plot(aux, main=NA)
  rug(panel$x)
  ##------------------------------------------------------------------
  ## faz o intervalo com o comprimento da banda
  arrows(panel$c0-0.5*panel$width, 0, panel$c0+0.5*panel$width, 0,
         length=0.1, code=3, angle=90, col=2)
  ##------------------------------------------------------------------
  ## faz uma seta que aponta para o valor da função
  y0 <- approx(aux$x, aux$y, xout=panel$c0)
  arrows(panel$c0, 0, panel$c0, y0$y, length=0.1, col=2)
  ##------------------------------------------------------------------
  ## desenha a função kernel acima do intevalo
  d <- density(panel$c0, width=panel$width, kernel=panel$kernel)
  lines(d$x, d$y/panel$n, col=2)
  ##------------------------------------------------------------------
  ## deve retornar a lista como última linha
  panel
  ##------------------------------------------------------------------
}

#-----------------------------------------------------------------------------
# alimenta os controladores, depois de rodar todos o painel gráfico está
# criado e então é só usar o mouse para alterar os componentes

# passar os argumentos que serão fixos, abre a janelinha
panel <- rp.control(x=rend, n=length(rend))

# controla a largura da banda
rp.slider(panel, width, 10, 1000, initval=150,
          showvalue=TRUE, action=density.panel)

# controla a coordenada do centro da banda
rp.slider(panel, c0, 3500, 6500, initval=5000,
          showvalue=TRUE, action=density.panel)

# controla o tipo de função kernel
rp.radiogroup(panel, kernel, 
              c("gaussian","epanechnikov","rectangular",
                "triangular","biweight","cosine","optcosine"),
              action=density.panel)

#-----------------------------------------------------------------------------
# obtendo valores de área abaixo da curva de densidade via aproximação da
# função densidade e integração numérica

dens <- density(rend, bw=150)
plot(dens)
rug(rend)
fdens <- approxfun(dens$x, dens$y) # cria uma função que é interpolação linear
i <- 4500; s <- 5500               # limites inferior e superior
integrate(fdens, i, s)             # faz uma integração numérica
polygon(c(i,dens$x[dens$x>i & dens$x<s],s),
        c(0,dens$y[dens$x>i & dens$x<s],0), col="gray90")
i <- 3600; s <- 6800               # limites inferior e superior
integrate(fdens, i, s)             # faz uma integração numérica
polygon(c(i,dens$x[dens$x>i & dens$x<s],s),
        c(0,dens$y[dens$x>i & dens$x<s],0), col="gray90")

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

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")
           )

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