Arquivo

Archive for the ‘matemática’ Category

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!

Fuzzy c-means, um exemplo!

Pois bem, por sugestão do Professor Julio da Motta Singer apresentarei um script para realizar uma análise de agrupamento no escopo da lógica fuzzy. Seja um experimento fatorial (a mesma natureza dos dados aqui disponibilizados) a seguinte implementação realiza a clusterização de uma dado fator marginalmente ao outro fator.

A heurística da lógica fuzzy tem um incremento interpretativo de, ao contrário do agrupamento k-means, que é dicotômico, associar a cada cluster uma “probabilidade” de associação. Todavia, no método arbitramos uma constante (m) de fuzzyficação. Assim o seguinte script realiza o agrupamento fuzzy (por intermédio da função cmeans() da biblioteca e1071) para vários valores m assim é possível observar o comportamento (decaimento) do erro associado a clusterização e também a disposição dos fatores em relação ao cluster.

Infelizmente (ou para alegria geral) não vou me atrever a discorrer da teoria envolvida nessa análise e um estudo específico deixo a cargo dos maiores interessados.

Todo esse mise en cene, me soa um pouco genérico demais, eu sei (–Me desculpem!). Mas acredito que o script possa ser útil de alguma forma.

Aproveito o ensejo para agradecer as sugestões do Prof Julio da Motta Singer ao trabalho e o incentivo a esse post. Também deixo um alô ao colegas que participaram, de uma forma ou de outra!

Vamos ao trabalho, primeiro instalando (para quem ainda não o fez) e carregando a biblioteca e o conjunto de dados!

## install.packages('e1071')
library(e1071) # carrega pacote com implementação 'c-means'

dados <- transform(read.table("http://dl.dropbox.com/u/38195533/fac_exemplo.txt", # arquivo
                              header = TRUE,
                              sep = '\t',
                              na.string = '.',
                              dec = ','),
                   fatorA = factor(fatorA),
                   fatorB = factor(fatorB))

str(dados)

ajuste <- lm(resp ~ fatorA + fatorB, data = dados)
anova(ajuste) # ANOVA

matriz <- tapply(dados$resp, list(dados$fatorB, dados$fatorA), mean)

A primeira etapa foi realizada. Não posso passar ou próximo trecho sem deixa-los a par do quem está por vir! O agrupamento que “tentamos” é com relação ao um padrão de comportamento do fator sob análise, em relação ao outro fator, assim constituímos ideótipos e assim procedemos. O modo como constituimos esses ideótipos é aplicado unicamente a natureza desse conjunto de dados, mas outros ideótipos, ao seu critério podem ser usados, fique a vontade!

Antes de passarmos agora a análise em si manipulamos os dados e constituímos os centróides:

FinWil <- apply(matriz, 2, mean) - mean(matriz) # Finlay & Wilkinson

altos <- which(FinWil >= 0) # ambientes acima da media
baixo <- which(FinWil < 0) # ambientes abaixo da media

medias <- cbind(apply(matriz[, altos], 1, mean),
                apply(matriz[, baixo], 1, mean))

##--------------------------------------------------------------------
## IDEOTIPOS
## I   - maximo em todos 
## II  - maximo nos altos minimo nos baixos
## III - minimo nos altos maximo nos baixos
## IV  - minimo em todos 
## ADICIONAIS:
## V   - media em todos
## VI  - maximo nos altos e media nos baixos
## VII  - media nos altos maximo nos baixos

I <- c(max(matriz[, altos]), max(matriz[, baixo]))
II <- c(max(matriz[, altos]), min(matriz[, baixo]))
III <- c(min(matriz[, altos]), max(matriz[, baixo]))
IV <- c(min(matriz[, altos]), min(matriz[, baixo]))
V <- c(mean(matriz[, altos]), mean(matriz[, baixo]))
VI <- c(max(matriz[, altos]), mean(matriz[, baixo]))
VII <- c(mean(matriz[, altos]), max(matriz[, baixo]))

centroides <- rbind(I, II, III, IV, V, VI, VII) # centróides

Resta agora realizar a análise em si, veja que o laço for() repete a análise para vários valores da constante m e de cada uma, armazena as “probabilidades” de associação a cada centróide e o erro do ajustamento!


ms <- seq(1.2, 10, by = .1)
lista <- array(NA, dim = c(nlevels(dados$fatorB), 
                           nrow(centroides), 
                           length(ms)))
erros <- matrix(NA, nrow = length(ms), ncol = 1)
dimnames(lista) <- list(levels(dados$TRAT), 
                        c('I', 'II', 'III', 'IV', 'V', 'VI', 'VII'), 
                        ms)

for(i in 1:length(ms))
  {
    j <- ms[i]
    ## função de cluster 'fuzzy c-means'
    grupamentos <- cmeans(medias, # matriz dos dados
                          centroides, # centróides dos clusters
                          verbose = FALSE, # não imprime histórico
                          dist = 'euclidean', # distância empregada
                          m = j, # constante de desfuzzyficação
                          method = 'cmeans') # método
    lista[,,i] <- grupamentos$membership
    erros[i,] <- grupamentos$withinerror
  }

Pronto… feito, quer dizer, quase! Realizamos portanto a análise para os vários valores de m resta agora “escolher” um deles e partir para o resto! Até a próxima…

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!

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.

Gerando seus próprios dados

Se você está prestes a ficar louco por causa dos seus dados que não ficaram bons ou porque você perdeu tudo! Esse post vem bem de encontro as suas expectativas.

Vamos mostrar agora como gerar seus próprios dados! Isso vai ser feito pelo uso das funções de geração de números aleatórios (pseudo aleatórios) do R, especialmente a rnorm().

É claro que é brincadeira que você vai usar os dados simulados como os dados verdadeiramente obtidos. Mas pretendemos com a demonstração da geração de dados mostrar uma boa funcionalidade do programa R.

E sem dúvida sempre que você abstrair um pouco para imaginar como gerar dados é um ótima forma de melhor compreender a mecânica envolvida no processo. Outra coisa é que a simulação é uma “bonita” forma de realizar testes, se você por acaso tiver uma idéia (ou hipótese) a testar mas ela seja meio difícil de se fazer experimentalmente a simulação é uma boa alternativa!

Para começar vamos apresentar a função rnorm(). É com ela que geramos os números pseudo-aleatórios. seus argumentos são n que é número de valores que serão gerados e depois vem os parâmetros da normal sob o qual os valores serão gerados mean e sd. Se você não especificar esses parâmetros serão gerado valores de uma distribuição normal padrão, media 0 e desvio padrão 1. O R tem função para geração de dados de outras distribuições, poisson, binomial, gamma, beta, etc. Procure a documentação!

Assim, primeiro temos que postular o que queremos simular. No nosso exemplo (CMR) vamos mostrar como gerar dados de um experimento em DIC, qualquer (qualquer mesmo…). Teremos então t tratamentos, r repetições e uma média geral mu, especificamos qual deva ser a variância entre tratamentos (sg.t) e uma variância do erro (sg.e).

Veja o CMR:

t <- 20 # número de tratamentos

r <- 10 # número de repetições

mu <- 100
sg.t <- 15
sg.e <- 9 # parâmetros (média geral, var de trat e var do erro)

a.dat <- matrix(rnorm(t, 0, sqrt(sg.t)),
                nrow = t,
                ncol = r) # iésimo efeito de tratamento

e.dat <- matrix(rnorm(t*r, 0, sqrt(sg.e)),
                nrow = t,
                ncol = r) # iésimo efeito aleatório

plan <- transform(data.frame(expand.grid(trat = 1:t, rep = 1:r), # dados
                            obs = c(round(mu + a.dat + e.dat))),
                 trat = factor(trat),
                 rep = factor(rep))

Esse conjunto de dados gerados é (digamos assim) “único”. A chance de você rodando novamente esse código obter exatamente os mesmos dados é bem improvável. Se você quiser repetir exatamente esses valores você pode especificar a semente que gera os dados, isso é feito pela função set.seed(). Faça o teste!

O modelo de DIC é: y_{ij} = \mu + t_i + \varepsilon_{ij} Portanto depois de especificar os parâmetros do nosso modelo criamos matrizes, com as dimensões do experimento que contêm simulações geradas com esses parâmetros.

Inspecione a matriz de nome a.dat, veja que a dimensão da matriz é t por r, mas são gerados apenas t simulações, assim todas células da mesma linha terão os mesmos valores, que serão os efeitos de tratamentos.

Por outro lado, a matriz e.dat tem todos os valores diferentes. Que são os desvios aleatórios. Como os efeitos de tratamentos e do erro são independentes somando-se as células dessas matrizes obtemos os dados.

É claro que esse é apenas uma simulação, é extremamente simples, mas ela serve como exemplo de como realizar uma simulação mais complexa, que envolva um outro delineamento diferente. E se você pretende usar a simulação para testar outras hipóteses seria interessante realizar mais de um conjunto de dados. E realizar alguns “testes” sobre os conjuntos gerados para verificar a consistência dos dados gerados. Mas vamos deixar isso para outro post!

Veja por exemplo que se você realizar uma análise sobre esse conjunto de dados as variâncias que você obterá não serão exatamente os parâmetros que você estipulou, isso ocorre pois como o exemplo tem uma dimensão pequena ocorrem alguns desvios devido ao processo de amostragem inerente.

Reproduzindo tabelas estatísticas

O objetivo desse “post” é dúbio! Vamos mostrar como reproduzir as tabelas estatísticas dos livros no R.

Mas para quê reproduzir as tabelas que estão nos livros? Primeiro vamos ver como fazer…

Primeiro pegue o seu livro de estatística experimental de cabeceira, e procure (no final, quase certamente) a tabela de F-Snedecor. Aquela mesma que agente usa para ver se um teste F< é significativo ou não! Se estivermos diante da mesma tabela teremos algo como (supondo que estejamos usando um nível de significância \alpha de 0,01, ou seja a tabela de 1%):

Resumo da Tabela dos Quantis da Distribuição F-Snedecor
gl_2 / gl_1 1 2 3 \cdots 10
1 4052,18 4999,50 5403,35 \cdots 6055,85
2 98,50 99,00 99,17 \cdots 99,40
3 34,12 30,82 29,46 \cdots 27,23
\vdots \vdots \vdots \vdots \ddots \vdots
100 6,90 4,82 3,98 \cdots 2,50

Para reproduzir essa tabela fazemos: criamos um vetor com os ditos numeradores (com os graus de liberdade) da fonte de variação a ser testada na análise de variância; depois criamos outro vetor com os graus de liberdade do denominador (o erro na maioria das vezes); com esses dois vetores usamos a função expand.grid() criando uma planilha com as margens da tabela de dupla entrada.

Isso feito (fácil né?). Agora obtemos os valores críticos (F_{tab}), que nada mais são que quantis da distribuição F-Snedecor. Com o vetor dos ditos F_{tab} concatenado a planilha que criamos anteriormente conseguimos “reproduzir” ipsis literis, as tabelas dos livros usando uma combinação das funções with() e tapply(), veja no CMR que reproduz a tabela de F!


alpha <- .01 
numerador <- 1:10 
denominador <- c(1:30, 35, 40, 45, 50, 100) 

margens <- expand.grid(numerador = numerador, # margens - tabela
                       denominador = denominador)

probs <- round(qf((1 - alpha),
                  margens$numerador,
                  margens$denominador), # quantil da dist. F-Snedecor
               dig = 2) # arredondamento

class <- transform(cbind(margens, probs), # planilha com as margens
                   numerador = factor(numerador), 
                   denominador = factor(denominador)) # transf. em fatores

tabela <- with(class,
               tapply(probs,
                      list(denominador, numerador),
                      sum)) # construção da tabela de dupla entrada

Você pode pensar “- Agora ficou fácil, nem preciso mais do livro, calculo os meus próprios F tabelados!”. Por trás desse “sentimento” de indepêndencia vamos agora fazer o mesmo para a tabela de \chi^2 (Veja o CMR). Não vamos nos alongar, mas acho que concordamos que isso é facilmente extrapolável para outras tabelas estatísticas do final dos livros.


gl <- c(1:30, 40, 50, 60, 120, 240, 480, 960)
nivel.prob <- sort(c(.005, .01, .025, .05, .1, 
                   seq(.25, .75, .25), .9, .95, .975, .990, .995),
                   decreasing = TRUE) 

margens <- expand.grid(gl = gl, alpha = nivel.prob) # margens - tabela de dupla entrada

chi2 <- qchisq((1 - margens$alpha),
               margens$gl) # quantil da distribuição de chi-quadrado

class <- transform(cbind(margens, chi2), # construção da planilha de valores
                   gl = factor(gl),
                   alpha = factor(alpha)) # transformação em fatores

tabela <- round(with(class,
               tapply(chi2,
                      list(gl, alpha),
                      sum)), dig = 3) # tabela de dupla entrada

Compare o que você obteve no R com a tabela do seu livro, vou reproduzir um “pedaço” dela aqui, caso você não esteja com seu livro:

Resumo da Tabela dos Quantis da Distribuição \chi^2
gl / \alpha 0,005 0,010 0,025 \cdots 0,995
1 7,879 6,635 5,024 \cdots <0,000
2 10,597 9,210 7,78 \cdots 0,010
3 12,838 11,345 9,348 \cdots 0,072
\vdots \vdots \vdots \vdots \ddots \vdots
960 1076,621 1064,867 1047,760 \cdots 850,891

Agora a lógica desse “post”: na verdade, reproduzir as tabelas foi só um subterfúgio. Até porque não tem coisa mais torpe que isso (as tabelas sempre continuarão nos livros). E principalmente, na verdade, elas, hoje em dia, perderam um pouco (para não dizer todo) o propósito! Antes de você ou seu professor mais ortodoxo de estatística me apedrejar, deixe eu terminar…!

Antigamente, obter esse valores de F_{tab} eram difíceis de serem obtidos! Isso envolve integração e cálculo “pesado”. Assim, é extremmamente louvável que algumas pessoas tiveram a intenção de fazer essas contas e construir essas tabelas. Mas hoje em dia, isso é facilmente obtido pelo computador (da mesma forma que fizemos!). E principalmente, abre-se um “leque” para nós, afinal não precisamos mais ficar “presos”, aos valores de \alpha das tabelas! Basta apenas que entendamos a relação com as probabilidades de erros…

Essa figura mostra as probabilidades de erro tipo I, tipo II e o poder do teste. Erro tipo I é o erro ao rejeitar H_0 quando, na realidade, H_0 é verdadeira. A probabilidade de cometer este erro do tipo I é designada por \alpha (nível de significância). O erro do tipo I equivale a concluir que o tratamentos diferem quando na verdade eles diferem.

A probabilidade de cometermos o erro tipo I é de \alpha e de cometermos o erro tipo II é \beta. Em um teste de hipótese ambas essas probabilidades devem ser mínimas. Mas em geral se “escolhe” pela diminuição do erro tipo I, que é feita manipulando-se \alpha.

grf-erros

A maioria do programas hoje em dia (incluindo o tema do Ridículas) apresentam além dos asteríscos da significância do teste F a estatística do p-valor que é nada mais nada menos que o complemento do da probabilidade do erro tipo I.

Categorias:dados, matemática Tags:, ,

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!