Arquivo

Posts Tagged ‘função’

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!

Anúncios

PCA de três dimensões

Olá a todos. O objetivo desse post é mostrar como fazer a análise de componentes principais em três dimensões. Mas esse post tem o destaque não ter sido escrito sozinho… Sem a participação da Msc Thalita (ou “Hermínia”, ninguém vai entender!) não teria post. Agradeçam a ela, mas reclamem comigo! Pois bem, Vamos usar como exemplo dados de um estudo com fatorial duplo (4 \times 4), em que foram tomadas quatro características (dados disponíveis para teste!).

Vou aproveitar e deixar aqui também um alô para o Dsc Guilherme pela ajuda nos detalhes do script, Valeu!

Um pouco de teoria:

Experimentos fatoriais (n \times n) caracterizam-se pela apresentação dos dados em tabelas de duas entradas (matriz), sendo que cada casela da tabela contém a resposta média de cada combinação dos fatores sob análise. Todas possibilidades desses experimentos estão fora do escopo desse post, por isso vamos um pouco além…

Acrescentando-se um outro fator, ou ainda outra característica (n \times n \times n, por exemplo), com 3 interações duplas e uma tripla, isso fica um pouco mais complexo. Neste caso, os dados são organizados em “CUBO” (arranjo de três entradas).

Nesta situação não é possível aplicar a análise de componentes principais usual. Uma alternativa é a utilização da decomposição PCAn, proposta por Tuckey em 1964.

Portanto nesse post pretendemos (afinal esse foi um trabalho em equipe) mostrar como realizar a decomposição de um arranjo de três entradas (CUBO). E mais que isso, mostrar uma função que realiza a escolha entre os vários modelos possíveis que variam conforme a dimensão dos dados.

Vamos ao código… mas calma, primeiro preparamos os dados!

Veja que usamos um pacote (PTAk, disponível no CRAN). Caso você não o tenha, proceda com install.packages(‘PTAk’)!

require(PTAk) # pacote analises pscicrometricas

## le e transforma os dados em fatores
dados <- transform(read.table("http://dl.dropbox.com/u/38195533/dados.txt")
                              na.string = '.',
                              header = TRUE,
                              dec = ',',
                              sep = '\t'),
                   elemento = factor(elemento),
                   residuo = factor(residuo),
                   metodo = factor(metodo))

CUBO.inicial <- with(dados,
                     tapply(conc, list(residuo, elemento, metodo),
                            mean, na.rm = TRUE))
## padronização dos dados
CUBO.standard <- apply(CUBO.inicial, 2, scale, center = FALSE)

## reconstrucao da planilha
estrutura <- expand.grid(elemento = factor(1:4),
                         residuo = factor(1:4),
                         metodo = factor(1:4))

plan <- data.frame(estrutura, conc = c(CUBO.standard)) ## dados

ajuste <- aov(conc ~ elemento * residuo * metodo - elemento:residuo:metodo,
              data = plan)
erros <- ajuste$residuals

## CUBO
Z <- tapply(erros, list(plan$residuo, plan$elemento, plan$metodo), mean) 

Nesse trecho são lidos os dados, na forma usual de uma planilha de dados (uma coluna para cada variável). Pela natureza dos dados, estes são padronizados (N(\mu = 0, \sigma^2 = 1)). Depois a panilha de dados é “reconstruída” e feito o ajuste do modelo linear respectivo (sem a interação tripla — objeto do nosso estudo!). Do resultado da análise extraímos os resíduos e construímos o CUBO.

Quase tudo pronto! O próximo trecho é a função que efetivamente faz os vários ajustes possíveis (isso é função da dimensão dos dados).

otimim.PCAn <- function(cubo) {
  dimensoes <- dim(Z) - 1
  diagnostico <- vector('list', length = prod(dimensoes - 1))
  ajustes <- vector('list', length = prod(dimensoes - 1))
  contador <- 0
  for(i in 2:dimensoes[ 1 ])
    {
      for(j in 2:dimensoes[ 2 ])
        {
          for(k in 2:dimensoes[ 3 ])
            {
              modelo <- c(i, j, k)
              tamanho <- sum(modelo)
              ## decomposicao Tucker -- PTAk
              tucker <- PCAn(Z, # array
                             dim = modelo, # dimensoes
                             test = 1E-12,
                             Maxiter = 1000,
                             smoothing = FALSE,
                             verbose = FALSE,
                             file = NULL)
              ## diagnostico
              porcentagem <- tucker[[ 3 ]]$pct # % explicada
              contador <- contador + 1
              ajustes[[ contador ]] <- tucker
              diagnostico[[ contador ]] <- c(modelo, tamanho, porcentagem)
            }
        }
    }
  ## diagnosticos dos ajustes
  tab.ajustes <- data.frame(do.call(rbind, diagnostico))
  tab.ajustes <- tab.ajustes[order(-tab.ajustes[, 4], -tab.ajustes[, 5]), ]
  por.total <- split(tab.ajustes, tab.ajustes[, 4])
  melhores <- vector('list', length = length(por.total))
  for(i in 1:length(melhores))
    {
      melhores[[i]] <- por.total[[i]][1, ]
    }
  ## melhores ajustes por dimensao
  tab.melhores <- do.call(rbind, melhores)
  tab.melhores <- tab.melhores[order(-tab.melhores[, 4]), ]
  ## melhor modelo
  melhor <- tab.melhores[(which((tab.melhores[,5]-
                          c(tab.melhores[,5][-1],0))>=5)[1]+1),]
  ## resposta -- melhor ajuste
  return(list(escolhido = melhor,
              resposta = ajustes[[ as.numeric(rownames(melhor)) ]]))
}

UFA!

Feita a função agora é só usa-lá (simples não?). Um uso possível dessa decomposição é representar e interpretar os dados graficamente em um biplot. Nessa condição de arranjos de três entradas são possíveis diferentes formas de gráficos biplots. Mas o objetivo do nosso post não é apresentá-los.

O trecho final aplica a função recém-construída (e retorna um controle de tempo de cada ajuste). Mostramos também como extrair as matrizes (A, B e C) e o arranjo G do melhor ajuste.

As matrizes de resposta e o arranjo G é que são úteis na construção dos gráficos. Aproveitem!

ajuste.PCA <- otimim.PCAn(Z)

A <- ajuste.PCA$resposta[[1]]$v
B <- ajuste.PCA$resposta[[2]]$v
C <- ajuste.PCA$resposta[[3]]$v

G <- ajuste.PCA$resposta[[3]]$coremat

summary(ajuste) # ANOVA do modelo
summary.PTAk(ajuste.PCA$resposta) # quanto cada componente explica

Voi lá…, Até mais pessoal!

Gráficos de funções e suas derivadas

Nessa dica vou apresentar o CMR (código mínimo reproduzível) para fazer gráfico de funções, obter derivadas simbólicas de funções, adicionar o gráfico das derivadas das funções com escala no eixo das abcissas da direita com as expressões na legenda do gráfico (ver figura abaixo).

Gráfico da função e sua primeira derivada.

 

A curve() é a função do R que faz gráficos no plano cartesiano de função paramétrica de uma variável. Diversas opções estão disponíveis para essa função, como os tipos de linhas, cores e demais parâmetros gráficos.

A D() é a função do R que retorna a expressão das derivadas (simbólicas) de uma função. O resultado desta e um objeto de classe language. Essa função funciona de maneira recursiva como você pode observar ao obter a derivada de segunda ordem. Infelizmente, o resultado da derivação não é simplificado matematicamente.

A sobreposição de gráficos é feita nas linhas [16,25,29] e a adição de uma escala para o eixo das abcissas foi colocada no lado direito do gráfico com o comando das linhas [18,28].

A função do.call() foi usada para tornar o procedimento mais automático no sentido de que com um objeto especificando a função obteremos todos os mesmos resultados. Saber automatizar os processos é importante quando tem que e avaliar diversas funções. Imagina ficar dando copia e cola o tempo todo? Vamos deixar isso para pessoas que usam programas (em geral pagos) que não oferecem os recursos que temos com o R. Tempo é dinheiro!

Com a função uniroot() obtemos as raízes da nossa função e representamos no gráfico com uma reta vertical.

Finalmente, confeccionamos uma figura para a publicação com rótulos em todos os eixos e legenda dentro do gráfico com a expressão de cada função.

#-----------------------------------------------------------------------------
# gráfico da função, um polinômio de ordem 3 com limites em x de -10 a 10

curve(1+1/2*x+1/3*x^2+1/4*x^3, -10, 10)

#-----------------------------------------------------------------------------
# derivada de primeira e segunda order dessa expressão com relação à x

D(expression(1+1/2*x+1/3*x^2+1/4*x^3), "x")
D(D(expression(1+1/2*x+1/3*x^2+1/4*x^3), "x") ,"x")

#-----------------------------------------------------------------------------
# sobre posição do gráfico da derivada primeira adicionando um novo eixo y

curve(1+1/2*x+1/3*x^2+1/4*x^3, -10, 10)
par(new=TRUE)
curve(1/2+1/3*(2*x)+1/4*(3*x^2), -10, 10, xlab="", ylab="", col=2, axes=FALSE)
axis(4, col=2)

#-----------------------------------------------------------------------------
# procedimento mais automático, atribuindo a expressão a um objeto

expr <- expression(1+1/2*x+1/3*x^2+1/4*x^3)
do.call(curve, list(parse(text=as.character(expr)), -10, 10, ylab=expr))
par(new=TRUE)
do.call(curve, list(D(expr, "x"), -10, 10,
                    xlab="", ylab="", col=2, axes=FALSE))
axis(4, col=2)
par(new=TRUE)
do.call(curve, list(D(D(expr, "x"), "x"), -10, 10,
                    xlab="", ylab="", col=3, axes=FALSE))

#-----------------------------------------------------------------------------
# adicionando a linha da raízes do polinômio

ro <- uniroot(function(x) 1+1/2*x+1/3*x^2+1/4*x^3, c(-10, 10))
str(ro)

curve(1+1/2*x+1/3*x^2+1/4*x^3, -10, 10)
abline(v=ro$root, col=1, h=0, lty=3)

#-----------------------------------------------------------------------------
# confeccionando uma figura para publicação, para salvar descomente

#png("f002.png", w=500, h=300);
par(mar=c(5.1,4.1,2.1,4.1))
expr <- expression(1+1/2*x+1/3*x^2+1/4*x^3)
do.call(curve, list(parse(text=as.character(expr)), -10, 10,
                    ylab=expression(f(x))))
par(new=TRUE)
do.call(curve, list(D(expr, "x"), -10, 10,
                    xlab="", ylab="", lty=2, axes=FALSE))
axis(4)
mtext(side=4, line=3, text=expression(df(x)/dx))
legend("topleft", bty="n", lty=1:2,
       legend=c(paste(expression(f(x)), "=", expr),
         paste(expression(df(x)/dx), "=", capture.output(D(expr, "x")))))
#dev.off()

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

Veja a documentação dessas funções para um melhor aproveitamento. Deixe suas sugestões/dúvidas nos comentários. Até a próxima Ridícula!