Apêndice BCálculos

Descrevemos brevemente dois pacotes computacionais que podem ser usados para todos os cálculos feitos no texto. Recomendamos que os estudantes se familiarizem com pelo menos um deles. A descrição do R é bastante completa, ao menos para os cálculos baseados no material deste texto, ao passo que é necessária outra referência para aprender o Minitab.

B.1 Usando o R

O R é um pacote estatístico gratuito que pode ser baixado e instalado no seu computador (veja http://cran.r-project.org/). Um manual gratuito também está disponível nesse site.

Depois de instalar o R no seu sistema, você pode iniciá-lo clicando no ícone correspondente (ou, em sistemas Unix, simplesmente digitando “R”). Aparece então uma janela, chamada R Console, que contém algum texto e um prompt “>”, depois do qual você digita os comandos. Os comandos são separados por novas linhas ou por “;”. A saída dos comandos também é exibida nessa janela, a menos que seja propositalmente direcionada para outro lugar. Para sair do R, digite q() depois do prompt. Para saber mais sobre qualquer coisa no R, um recurso conveniente é usar o Help da barra de menus, no alto da janela do R. Outra opção é digitar ?name depois do prompt (e pressionar Enter) para exibir informações sobre name; por exemplo, ?q abre uma página com informações sobre o comando de encerramento q.

Operações e Funções Básicas

Um comando básico avalia uma expressão, como

> 2+3
[1] 5

que soma 2 e 3 e produz a resposta 5. Outra possibilidade é atribuir o valor da expressão a uma variável, como em

> a <- 2

em que <- (sinal de menor seguido de sinal de menos) atribui o valor 2 a uma variável chamada a. Também se pode usar = na atribuição, como em a = 2, mas usaremos <-. Podemos então verificar essa atribuição simplesmente digitando a e apertando Enter, o que faz o valor de a ser impresso.

> a
[1] 2

Observe que o R diferencia maiúsculas de minúsculas, de modo que A seria uma variável diferente de a. Há algumas restrições na escolha de nomes para variáveis e vetores, mas você não errará se sempre começar o nome com uma letra.

Podemos atribuir os valores de um vetor usando a função de concatenação c(), como em

> b <- c(1,1,1,3,4,5)
> b
[1] 1 1 1 3 4 5

o que cria um vetor chamado b com seis valores. Podemos acessar a \(i\)-ésima entrada de um vetor b referindo-nos a ela como b[i]. Por exemplo,

> b[3]
[1] 1

imprime a terceira entrada de b, isto é, 1. Outra possibilidade é usar o comando scan para a entrada de dados. Por exemplo,

> b <- scan()
1:  1 1 1 3 4 5
7:
Read 6 items
> b
[1] 1 1 1 3 4 5

faz a mesma atribuição. Observe que, com o comando scan, simplesmente digitamos os dados e encerramos a entrada com uma linha em branco. Também podemos usar scan para ler dados de um arquivo; para isso, remetemos o leitor a ?scan.

Às vezes queremos vetores cujas entradas seguem algum padrão. Muitas vezes podemos usar a função rep para isso. Por exemplo, x <- rep(1,20) cria um vetor de 20 uns. Padrões mais complicados podem ser obtidos; para isso, remetemos o leitor a ?rep.

A aritmética básica pode ser feita com variáveis e vetores usando + (adição), - (subtração), * (multiplicação), / (divisão) e ^ (exponenciação). Essas operações são feitas componente a componente. Por exemplo, poderíamos multiplicar cada componente de b por ele mesmo com

> b*b
[1] 1 1 1 9 16 25

ou multiplicar cada elemento de b por 2, como em

> 2*b
[1] 2 2 2 6 8 10

o que realiza essa operação.

Há várias funções disponíveis no R, como abs(x) (calcula o valor absoluto de x), log(x) (calcula o logaritmo natural de x), exp(x) (calcula \(e\) elevado à potência x), sin(x), cos(x), tan(x) (que calculam as funções trigonométricas), sqrt(x) (que calcula a raiz quadrada de x), ceiling(x) e floor(x) (que calculam o teto e o piso de x). Quando uma dessas funções é aplicada a um vetor x, ela devolve um vetor de mesmo comprimento, com a função aplicada a cada elemento do vetor original. Há inúmeras funções especiais disponíveis no R, mas duas importantes são gamma(x), que devolve a função gama¹ aplicada a x, e lgamma(x), que devolve o logaritmo natural da função gama.

Há também funções que devolvem um único valor quando aplicadas a um vetor. Por exemplo, min(x) e max(x) devolvem, respectivamente, o menor e o maior elemento de x; length(x) dá o número de elementos de x; e sum(x) dá a soma dos valores de x.

O R também opera com as quantidades lógicas TRUE (ou T, para verdadeiro) e FALSE (ou F, para falso). Valores lógicos são gerados por condições que são verdadeiras ou falsas. Por exemplo,

> a <- c(-3,4,2,-1,-5)
> b <- a>0
> b
[1] FALSE TRUE TRUE FALSE FALSE

compara cada elemento do vetor a com 0, devolvendo TRUE quando ele é maior que 0 e FALSE caso contrário, e esses valores lógicos são guardados no vetor b. Podem ser usados os seguintes operadores lógicos: <, <=, >=, >, == (igualdade), != (desigualdade), além de & (conjunção), | (disjunção) e ! (negação). Por exemplo, se criarmos um vetor lógico c assim:

> c <- c(T,T,T,T,T)
> b&c
[1] FALSE TRUE TRUE FALSE FALSE
> b|c
[1] TRUE TRUE TRUE TRUE TRUE

então um elemento de b&c é TRUE quando os dois elementos correspondentes de b e c são TRUE, enquanto um elemento de b|c é TRUE quando pelo menos um dos elementos correspondentes de b e c é TRUE.

Às vezes podemos ter variáveis que assumem valores de caractere. Embora sempre seja possível codificar esses valores como números, não há necessidade disso, pois o R também lida com variáveis de valores de caractere. Por exemplo, os comandos

> A <- c('a','b')
> A
[1] "a" "b"

criam um vetor de caracteres A, contendo dois valores, a e b, e depois imprimem esse vetor. Observe que pusemos os valores de caractere entre aspas simples ao fazer a atribuição.

Às vezes faltam valores nos dados, e eles são então listados como NA (not available, não disponível). Operações com valores faltantes criam valores faltantes. Além disso, uma operação impossível, como 0/0, produz NaN (not a number, não é um número).

Vários objetos podem ser criados durante uma sessão do R. Para ver os que foram criados até o momento na sua sessão, use o comando ls(). Você pode remover qualquer objeto da sua área de trabalho usando o comando rm. Por exemplo, rm(x) remove o vetor x.

Funções de Probabilidade

O R tem várias funções embutidas para calcular a fda, a fda inversa, a função densidade² ou de probabilidade e para gerar amostras aleatórias das distribuições comuns que encontramos em probabilidade e estatística. Elas se distinguem por um prefixo e pelo nome da distribuição base. Alguns dos nomes de distribuições são dados na tabela a seguir.

Como sempre, é preciso tomar cuidado com a distribuição gama. O caminho mais seguro é incluir outro argumento junto com a distribuição para indicar se lambda é um parâmetro de taxa (a densidade é \((\Gamma(\alpha))^{-1}\lambda^{\alpha}x^{\alpha-1}e^{-\lambda x}\)) ou um parâmetro de escala (a densidade é \((\Gamma(\alpha))^{-1}\lambda^{-\alpha}x^{\alpha-1}e^{-x/\lambda}\)). Assim, gamma(\(\cdot\),alpha,rate=lambda) indica que lambda é um parâmetro de taxa, e gamma(\(\cdot\),alpha,scale=lambda) indica que ele é um parâmetro de escala.

O argumento indicado por \(\cdot\) é especificado de acordo com a finalidade do comando que usa o nome da distribuição. Para obter a fda de uma distribuição, ponha p antes do nome; então \(\cdot\) é o valor em que se quer calcular a fda. Para obter a fda inversa de uma distribuição, ponha q antes do nome; então \(\cdot\) é o valor em que se quer calcular a fda inversa. Para obter a função densidade ou de probabilidade, ponha d antes do nome; então \(\cdot\) é o valor em que se quer calcular a densidade ou a função de probabilidade³. Para obter amostras aleatórias, ponha r antes do nome; então \(\cdot\) é o tamanho da amostra⁴ aleatória que se quer gerar.

Por exemplo,

> x <- rnorm(4,1,2)
> x
[1] -0.2462307 2.7992913 4.7541085 3.3169241

gera uma amostra de tamanho 4 da distribuição \(N(1, 2^2)\) e a atribui ao vetor x. O comando

> dnorm(3.2,2,.5)
[1] 0.04478906

calcula a densidade \(N(2; 0{,}25)\) em 3,2, enquanto

> pnorm(3.2,2,.5)
[1] 0.9918025

calcula a fda \(N(2; 0{,}25)\) em 3,2, e

> qnorm(.025,2,.5)
[1] 1.020018

dá o quantil 0,025 da distribuição \(N(2; 0{,}25)\).

Se temos dados guardados num vetor x, podemos amostrar valores de x, com ou sem reposição, usando a função sample. Por exemplo, sample(x,n,T) gera uma amostra de tamanho n de x com reposição, enquanto sample(x,n,F) gera uma amostra de tamanho n de x sem reposição (observe que, neste último caso, n não pode ser maior que length(x)).

Às vezes é conveniente poder repetir uma simulação de modo que os mesmos valores aleatórios sejam gerados. Para isso, pode-se usar o comando set.seed. Por exemplo, set.seed(12345) fixa a semente em 12345.

Tabulação de Dados

O comando table está disponível para tabular dados. Por exemplo, table(x) devolve uma tabela com a lista dos valores distintos encontrados em x e suas frequências de ocorrência em x. Essa tabela pode ser atribuída a uma variável com

y <- table(x)

para análise posterior (veja o trecho O Teste Qui-Quadrado, mais adiante).

Se x e y são vetores de mesmo comprimento, então table(x,y) produz uma tabulação cruzada, isto é, conta o número de vezes em que se obtém cada valor possível de \((x, y)\), em que \(x\) pode ser qualquer um dos valores assumidos em x e \(y\) pode ser qualquer um dos valores assumidos em y.

Gráficos de Dados

O R tem vários comandos disponíveis para fazer gráficos de dados. Por exemplo, suponha que temos uma amostra de tamanho \(n\) guardada no vetor x.

O comando hist(x) fornece um histograma de frequências dos dados, com os pontos de corte escolhidos automaticamente pelo R. Podemos acrescentar argumentos opcionais a hist. Seguem alguns dos argumentos disponíveis.

Por exemplo, hist(x,breaks=c(-10,-5,-2,0,2,5,10),freq=F) faz o gráfico de um histograma de densidade com pontos de corte \(-10, -5, -2, 0, 2, 5, 10\), em que tivemos o cuidado de garantir que min(x) \(> -10\) e max(x) \(< 10\).

Se y é outro vetor de mesmo comprimento que x, podemos produzir um diagrama de dispersão⁵ de y contra x com o comando plot(x,y). O comando plot(x,y,type="l") fornece um diagrama de dispersão de y contra x, mas agora com os pontos ligados por linhas. O comando plot(x) faz o gráfico dos valores de x contra seus índices. O comando plot(ecdf(x)) faz o gráfico da fda empírica dos dados de x.

Um boxplot dos dados de x é obtido com o comando boxplot(x). Boxplots lado a lado dos dados de x, y, z etc. podem ser obtidos com boxplot(x,y,z).

Um gráfico de probabilidade normal⁶ dos valores de x pode ser obtido com o comando qqnorm(x).

Um gráfico de barras pode ser obtido com o comando barplot. Por exemplo,

> h <- c(1,2,3)
> barplot(h)

produz um gráfico de barras com 3 barras de alturas 1, 2 e 3.

Há muitos outros aspectos dos gráficos no R que dão ao usuário um controle considerável sobre a aparência dos gráficos. Remetemos o leitor ao manual para uma discussão mais ampla sobre eles.

Inferência Estatística

O R tem uma abordagem poderosa para ajustar modelos e fazer inferência sobre eles. Os modelos são especificados pelo símbolo ~. Não discutimos isso por completo aqui; apenas indicamos como usá-lo para tratar os modelos de regressão linear simples e múltipla (em que a resposta e os preditores são todos quantitativos), os modelos com um e com dois fatores (em que a resposta é quantitativa, mas os preditores são categóricos) e o modelo de regressão logística⁷ (em que a resposta é categórica, mas os preditores são quantitativos). Suponha, então, que temos um vetor y com os valores da resposta.

Estatísticas Básicas

A função mean(y) devolve a média dos valores de y, var(y) devolve a variância⁸ amostral dos valores de y, e sd(y) dá o desvio padrão⁹ amostral. O comando median(y) devolve a mediana de y, enquanto quantile(y,p) devolve os quantis amostrais especificados no vetor de probabilidades p. Por exemplo, quantile(y,c(.25,.5,.75)) devolve a mediana e o primeiro e o terceiro quartis. A função sort(y) devolve um vetor com os valores de y ordenados do menor para o maior, e rank(y) dá os postos dos valores de y.

O Teste \(t\)

Para os dados em y, podemos usar o comando

> t.test(y,mu=1,alternative="two.sided",conf.level=.95)

para fazer um teste \(t\). Ele calcula o valor-p do teste de \(H_0: \mu = 1\) e forma um intervalo de confiança¹⁰ 0,95 para \(\mu\).

O Teste Qui-Quadrado

Suponha que y contenha um vetor de contagens para k caselas e que prob contenha as probabilidades hipotéticas dessas caselas. Então o comando

> chisq.test(y,p=prob)

faz o teste qui-quadrado para avaliar essa hipótese. Observe que y também poderia corresponder a uma tabela unidimensional.

Se x e y são dois vetores de mesmo comprimento, então chisq.test(x,y) faz um teste qui-quadrado de independência na tabela formada pela tabulação cruzada das entradas de x e y. Se primeiro criarmos essa tabulação cruzada na tabela t usando a função table, então chisq.test(t) faz esse teste.

Regressão Linear Simples

Suponha que temos um único preditor¹¹, com valores no vetor x. O modelo de regressão linear simples¹² \(E(y \mid x) = \beta_1 + \beta_2 x\) é então especificado no R por y~x. Chamamos y~x de fórmula do modelo e a lemos como “y é modelado como um modelo linear envolvendo x”. Para fazer o ajuste (o que fizemos aqui para um conjunto específico de dados), usamos o comando de ajuste de modelos lineares lm, como segue. O comando

> regexamp <- lm(y~x)

faz os cálculos de ajuste e de inferência sobre esse modelo e atribui o resultado a uma estrutura chamada regexamp. Qualquer outro nome válido poderia ter sido usado para essa estrutura. Podemos agora usar várias funções do R para extrair vários itens de interesse. Por exemplo,

> summary(regexamp)
Call:
lm(formula = y~x)
Residuals:
  Min      1Q    Median 3Q     Max
-4.2211 -2.1163 0.3248 1.7255 4.3323
Coefficients:
 Estimate          Std.  Error  t value Pr(>|t|)
(Intercept) 6.5228   1.2176     5.357   4.31e-05 ***
    x       1.7531   0.1016    17.248   1.22e-12 ***
--
Signif.  codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’
0.1 ‘ ’ 1
Residual standard error:  2.621 on 18 degrees of freedom
Multiple R-squared:  0.9429, Adjusted R-squared:  0.9398

F-statistic:  297.5 on 1 and 18 DF, p-value:  1.219e-12

usa a função summary para nos dar todas as informações de que precisamos. Por exemplo, a reta ajustada é dada por \(6{,}5228 + 1{,}7531x\). O teste de \(H_0: \beta_2 = 0\) tem valor-p de \(1{,}22 \times 10^{-12}\), de modo que temos forte evidência contra \(H_0\). Além disso, o \(R^2\) é dado por 94,29%. Itens individuais podem ser acessados por várias funções do R; para isso, remetemos o leitor a ?lm.

Regressão Linear Múltipla

Se temos dois preditores quantitativos nos vetores x1 e x2, podemos proceder exatamente como na regressão linear simples para ajustar o modelo de regressão linear \(E(y \mid x_1, x_2)\) \({}= \beta_1 + \beta_2 x_1\) \({}+ \beta_2 x_2\).¹³ Por exemplo, os comandos

> regex <- lm(y~x1+x2)
> summary(regex)

ajustam o modelo linear acima, atribuem os resultados à estrutura regex, e então a função summary imprime (omitido aqui) todas as quantidades relevantes. Lemos y~x1+x2 como “y é modelado como um modelo linear envolvendo x1 e x2”. Em particular, obtêm-se a estatística \(F\) e o valor-p associado a ela para o teste de \(H_0: \beta_2 = \beta_3 = 0\).

Isso se generaliza imediatamente para modelos de regressão linear com \(k\) preditores quantitativos \(x_1, \ldots, x_k\). Além disso, suponha que queremos testar se o modelo envolve apenas \(x_1, \ldots, x_l\), para \(l < k\). Usamos lm para ajustar o modelo com todos os \(k\) preditores e atribuímos isso a regex, e também usamos lm para ajustar o modelo que envolve só \(l\) preditores e atribuímos isso a regex1. Então o comando anova(regex,regex1) fornece a estatística \(F\) e seu valor-p para o teste de \(H_0: \beta_{l+1} = \cdots = \beta_k = 0\).

ANOVA com Um e com Dois Fatores

Suponha agora que \(A\) denota um preditor categórico que assume dois níveis, \(a_1\) e \(a_2\). Observe que os valores de \(A\) podem ser de caractere em vez de numéricos; por exemplo, x é um vetor de caracteres contendo os valores a1 e a2, usados para indicar em que nível o valor correspondente de y foi observado. Em qualquer caso, precisamos transformar isso num fator A, com o comando

> A <- factor(x)

para que A possa ser usado na análise. Então o comando

> aov(y~A)

produz a tabela de ANOVA com um fator. Naturalmente, aov também trata fatores com mais de dois níveis. Para obter as médias das caselas, use o comando tapply(y,A,mean).

Suponha que exista um segundo fator \(B\), com 5 níveis \(b_1, \ldots, b_5\). Se esse é o fator B no R, então o comando

> aov(y~A+B+A:B)

produz a ANOVA com dois fatores para testar interações entre os fatores A e B. Para obter as médias das caselas, use o comando tapply(y,list(A,B),mean). O comando aov(y~A+B) produz a tabela de ANOVA supondo que não há interações.

Regressão Logística

Suponha que temos dados binários guardados no vetor y e que x contém os valores correspondentes de um preditor quantitativo. Então podemos usar o comando de modelo linear generalizado glm para ajustar o modelo de regressão logística \(P(Y = 1 \mid x)\) \({}= \exp\{\beta_1\) \({}+ \beta_2 x\}/(1 + \exp\{\beta_1 + \beta_2 x\})\). Os comandos

> logreg <- glm(y~x,family=binomial)
> summary(logreg)

ajustam o modelo de regressão logística, atribuem os resultados a logreg, e então o comando summary imprime esse material. Isso nos dá as estimativas dos \(\beta_i\), seus erros padrão e os valores-p para testar que \(\beta_i = 0\).

Comandos de Controle e Programas em R

Um comando de controle básico tem a forma if (expr1) expr2 else expr3, em que expr1 assume um valor lógico, expr2 é executado se expr1 é T, e expr3 é executado se expr1 é F. Por exemplo, se x é uma variável que assume o valor \(-2\), então

> if (x<0) {y <- -1} else {y <- 1}

faz com que seja atribuído a y o valor \(-1\). Observe que a parte else do comando pode ser omitida.

O comando for (name in expr) expr2 executa expr2 para cada valor de name em expr1.¹⁴ Por exemplo,

> for (i in 1:10) print(i)

imprime o valor da variável i à medida que são atribuídos a i, sucessivamente, os valores de \(\{1, 2, \ldots, 10\}\). Observe que m:n é, no R, uma abreviação da sequência \((m, m + 1, \ldots, n)\). Outro exemplo:

> for (i in 1:20) y[i] <- 2^i

cria um vetor y com 20 entradas, em que o \(i\)-ésimo elemento de y é igual a \(2^i\).

O break encerra um laço, possivelmente com base em alguma condição que se verifique, enquanto next interrompe o processamento da iteração atual e avança o índice do laço. Tanto break quanto next valem apenas para o mais interno de laços aninhados.

Comandos no R podem ser agrupados colocando-os entre chaves, {expr1; expr2; ...}. Os comandos dentro das chaves são executados como uma unidade. Por exemplo,

> for (i in 1:20) {print(i); y[i] <- 2^i}; print(y[i])}

faz com que i seja impresso, y[i] seja atribuído e y[i] seja impresso, tudo dentro de um laço for.¹⁵

Muitas vezes, quando um cálculo é complicado, como um que envolve laços, é melhor pôr todos os comandos do R num único arquivo e executar o arquivo em modo de lote. Por exemplo, suponha que você tem um arquivo prog.R contendo código R. Então o comando source("pathname/prog.R") faz com que todos os comandos do arquivo sejam executados.

Muitas vezes é conveniente pôr comentários nos programas em R para explicar o que as linhas de código fazem. Uma linha de comentário é precedida por # e, claro, não é executada.

Funções Definidas pelo Usuário

O R também permite funções definidas pelo usuário. A sintaxe da definição de uma função é a seguinte.

Por exemplo, a função a seguir calcula o coeficiente de variação amostral dos dados x.

coef_var <- function(x) {
    result <- sd(x)/mean(x);
    return(result);
}

Então, se depois quisermos calcular o coeficiente de variação dos dados y, basta digitar coef_var(y).

Arranjos e Listas

Um vetor de comprimento \(m\) também pode ser visto como um arranjo unidimensional de comprimento \(m\). O R pode lidar com arranjos multidimensionais, por exemplo, arranjos \(m \times n\), \(m \times n \times p\) etc. Se a é um arranjo tridimensional, então a[i,j,k] se refere à entrada na posição \((i, j, k)\) do arranjo. Há várias operações que podem ser feitas com arranjos; para elas, remetemos o leitor ao manual. Mais adiante neste manual, discutiremos o caso especial dos arranjos bidimensionais, também conhecidos como matrizes. Por ora, pensamos nos arranjos apenas como objetos em que guardamos dados.

Uma estrutura de dados muito geral no R é dada por uma lista. Uma lista é parecida com um arranjo, com várias diferenças importantes.

1. Qualquer entrada de um arranjo é referida pelo seu índice. Já qualquer entrada de uma lista pode ser referida por um nome de caracteres. Por exemplo, os coeficientes da regressão ajustada são referidos por regex$coefficients depois do ajuste do modelo linear regex <- lm(y ~ x1 + x2). O cifrão ($) é o operador de referência a entradas, isto é, varname$entname indica a entrada “entname” da lista “varname”.
2. Enquanto um arranjo guarda dados de um só tipo, uma lista pode guardar quaisquer objetos do R. Por exemplo, a entrada coefficients de um objeto de regressão linear é um vetor numérico, e a entrada model é uma lista.
3. Os operadores de referência são diferentes: arr[i] se refere à \(i\)-ésima entrada do arranjo arr, e lst[[i]] se refere à \(i\)-ésima entrada da lista lst. Observe que \(i\) pode ser o nome da entrada, isto é, lst$entname e lst[['entname']] se referem aos mesmos dados.

Exemplos

Consideramos agora alguns exemplos relativos a seções ou exemplos específicos do texto principal. Para executar qualquer um destes códigos, é preciso primeiro definir as funções. Para isso, carregue o código com o comando source. Depois, é preciso especificar os argumentos das funções. Observe que as linhas das listagens podem estar quebradas de modo pouco natural e continuar na linha seguinte.

#   Example B.1.1
#   function name: bootstrap_median
#   parameters:
#       m   resample size
#       x   original data
#   return value:
#       a vector of resampled medians
#   description: resamples and stores its median
  bootstrap_median <- function(m,x) {
      n <- length(x);
      result <- rep(0,m);
      for(i in 1:m) result[i] <- median(sample(x,n,T));
   return(result);
  }

#   Example B.1.2
#   function name: post_normal
#   parameters:
#       m       sample size
#       x       data
#       alpha0  shape parameter for 1/sigma^2
#       beta0   rate parameter for 1/sigma^2
#       mu0     location parameter for mu
#       tau0square  variance ratio parameter for mu
#   returned values:
#       mu      sampled mu
#       sigmasq sampled sigmasquare
#   description: samples from the posterior distribution
#   in Example 7.3.1
#
post_normal<-function(m,x,alpha0,beta0,mu0,tau0square){
# set the length of the data
n <- length(x);
# the shape and rate parameters of the posterior dist.
# alpha_x = first parameter of the gamma dist.
#         = (alpha0 + n/2)
alpha_x <- alpha0 + n/2
# beta_x  = the rate parameter of the gamma dist.
beta_x  <- beta0 + (n-1)/2 * var(x) + n*(mean(x)-mu0)**2/
           2/(1+n*tau0square);
# mu_x    = the mean parameter of the normal dist.
mu_x    <- (mu0/tau0square+n*mean(x))/(n+1/tau0square);
# tausq_x = the variance ratio parameter of the normal
#           distribution
tausq_x <- 1/(n+1/tau0square);
# initialize the result
result  <- list();
result$sigmasq <- 1/rgamma(m,alpha_x,rate=beta_x);
result$mu      <- rnorm(m,mu_x,sqrt(tausq_x *
                      result$sigmasq));
return(result);
}

# Example B.1.3
# set the data
x <- c(11.6714, 1.8957, 2.1228, 2.1286, 1.0751, 8.1631,
  1.8236, 4.0362, 6.8513, 7.6461, 1.9020, 7.4899,
  4.9233, 8.3223, 7.9486);
post <- post_normal(10**4,x,2,1,4,2);
# compute the coefficient of variation
psi <- sqrt(post$sigmasq)/post$mu;
psi_hat <- mean(psi <= .5);
psi_se  <- sqrt(psi_hat * (1-psi_hat))/sqrt(length(psi));
# the interval
cq <- 3
cat("The three times s.e. interval is ",
  "[",psi_hat-cq*psi_se, ", ", psi_hat+cq*psi_se,"]\n");

# Example B.1.4
#   function name: gibbs_normal
#   parameters
#     m         the size of posterior sample
#     x         data
#     alpha0    shape parameter for 1/sigma^2
#     beta0     rate parameter for 1/sigma^2
#     lambda    degree of freedom of Student's t-dist.
#     mu0       location parameter for mu
#     tau0sq    scale parameter for mu
#     burnin    size of burn in. the default value is 0.
#
#   returnrd values
#       mu       sampled mu's
#       sigmasq  sampled sigma^2's
#   description: samples from the posterior in Ex. 7.3.2
#
gibbs_normal <- function(m,x,alpha0,beta0,lambda,mu0,
                         tau0sq,burnin=0) {
  # initialize the result
  result <- list();
  result$sigmasq <- result$mu <- rep(0,m);
  # set the initial parameter
  mu      <- mean(x);
  sigmasq <- var(x);
  n       <- length(x);

  # set parameters
  alpha_x <- n/2 + alpha0 + 1/2;
  # loop
  for(i in (1-burnin):m) {
   # update v_i's
   v       <- rgamma(n,(lambda+1)/2,rate=((x-mu)**2/
                          sigmasq/lambda+1)/2);
   # update sigma-square
   beta_x <-(sum(v*(x-mu)**2)/lambda+(mu-mu0)**2/
                               tau0sq)/2+beta0;
   sigmasq<- 1/rgamma(1,alpha_x,rate=beta_x);
   # update mu
   r       <- 1/(sum(v)/lambda+1/tau0sq);
   mu      <- rnorm(1,r*(sum(v*x)/lambda+mu0/tau0sq),
                       sqrt(r*sigmasq));
   # burnin check
   if(i < 1) next;
   result$mu[i]      <- mu;
   result$sigmasq[i] <- sigmasq;
  }
  result$psi <- sqrt(result$sigmasq)/result$mu;
  return(result);
}

#   Example B.1.5
#   function name: batching
#   parameters:
#       k   size of each batch
#       x   data
#   return value:
#       an array of the averages of each batch
#   description: this function separates the data x into
#   floor(length(x)/k) batches and returns the array of
#   the averages of each batch
batching <- function(k,x) {
    m <- floor(length(x)/k);
    result <- rep(0,m);
    for(i in 1:m) result[i] <- mean(x[(i-1)*k+(1:k)]);
    return(result);
}

#   Example B.1.6
#   function name: discrepancy
#   parameters:
#       m   resample size
#       n   size of data
#   return value:
#       an array of m discrepancies
#   description: this function generates m discrepancies
#   when the data size is n
discrepancy <- function(m,n) {
    result <- rep(0,m);
    for(i in 1:m) {
        x    <- rnorm(n);
        xbar <- mean(x);
        r    <- (x-xbar)/sqrt((sum((x-xbar)**2)));
        result[i] <- -sum(log(r**2));
    }
    return(result/n);
}

#   Example B.1.7
#   function name: ddirichlet
#   parameters:
#       n      sample size
#       alpha  vector(alpha1,...,alphak)
#   return value:
#       a (n x k) matrix. rows are i.i.d. samples
#   description: this function generates n random samples
#   from Dirichlet(alpha1,...,alphak) distribution
ddirichlet <- function(n,alpha) {
    k <- length(alpha);
    result <- matrix(0,n,k);
    for(i in 1:k) result[,i] <- rgamma(n,alpha[i]);
    for(i in 1:n) result[i,] <- result[i,] /
                                     sum(result[i,]);
    return(result);
}

Matrizes

Uma matriz pode ser vista como uma coleção de valores de dados com dois índices ou como um arranjo retangular de dados. Assim, se a é uma matriz, então a[i,j] é o elemento \((i, j)\) de a. Observe que a[i,] se refere à \(i\)-ésima linha de a e a[,j] se refere à \(j\)-ésima coluna de a. Se uma matriz tem \(m\) linhas e \(n\) colunas, ela é uma matriz \(m \times n\), e \(m\) e \(n\) são chamados de dimensões da matriz.

Talvez a maneira mais simples de criar matrizes seja com os comandos cbind e rbind. Por exemplo,

> x<-c(1,2,3)
> y<-c(4,5,6)
> a<-cbind(x,y)
> a
 x y
[1,] 1 4
[2,] 2 5
[3,] 3 6

cria os vetores x e y, e o comando cbind toma x como primeira coluna e y como segunda coluna da matriz \(3 \times 2\) recém-criada a. Observe que, na impressão de a, as colunas continuam rotuladas x e y, embora ainda possamos nos referir a elas como a[,1] e a[,2]. Podemos remover esses nomes de colunas com o comando colnames(a)<-NULL. Do mesmo modo, o comando rbind trata os vetores dados como argumentos como as linhas de uma matriz. Para determinar o número de linhas e de colunas de uma matriz a, podemos usar os comandos nrow(a) e ncol(a). Também podemos criar uma matriz diagonal usando o comando diag. Se x é um vetor \(n\)-dimensional, então diag(x) é uma matriz \(n \times n\) com as entradas de x na diagonal e 0 nas demais posições. Se a é uma matriz \(m \times n\), então diag(a) é o vetor com as entradas tiradas da diagonal principal de a. Para criar uma matriz identidade \(n \times n\), use diag(n).

Há várias operações que podem ser feitas com matrizes. Se as matrizes a e b são \(m \times n\), então a+b é a matriz \(m \times n\) formada pela soma das matrizes componente a componente. A transposta de a é a matriz \(n \times m\) t(a), cuja \(i\)-ésima linha é igual à \(i\)-ésima coluna de a. Se c é um número, então c*a é a matriz \(m \times n\) formada pela multiplicação de cada elemento de a por c. Se a é \(m \times n\) e b é \(n \times p\), então a%*%b é a matriz \(m \times p\) produto (Apêndice A.4) de a e b. Um vetor numérico é tratado como vetor coluna na multiplicação de matrizes. Observe que a*b também está definido quando a e b têm a mesma dimensão, mas esse é o produto componente a componente das duas matrizes, bem diferente do produto de matrizes.

Se a é uma matriz \(m \times m\), a inversa de a é obtida como solve(a). O comando solve devolve um erro se a matriz não tiver inversa. Se a é uma matriz quadrada, então det(a) calcula o determinante de a.

Consideramos agora uma aplicação importante.

Pacotes

Há muitos pacotes escritos para estender a capacidade do R básico. É muito provável que, se você tiver uma necessidade de análise de dados que não possa ser atendida pelo R, encontre um pacote disponível gratuitamente para acrescentar. Remetemos o leitor a ?install.packages e ?library para saber mais sobre isso.

B.2 Usando o Minitab

Todos os cálculos encontrados neste texto foram feitos com o Minitab. Esse pacote de software estatístico é muito fácil de aprender e de usar. Outros pacotes, como o SAS ou o R (ver a Seção B.1), também poderiam ser usados com essa finalidade.

A maior parte dos cálculos foi feita usando o Minitab como uma calculadora, isto é, os dados eram digitados e então se recorria a vários comandos do Minitab para obter as quantidades desejadas. Esses cálculos não exigem nenhuma programação.

Houve alguns cálculos, porém, que envolveram um pouco de programação. Em geral, tratava-se de um cálculo em que muitas operações tinham de ser feitas repetidas vezes, e por isso convinha usar um laço. Em cada um desses casos, registramos aqui o código do Minitab que usamos nesses cálculos. Como mostram os exemplos a seguir, esses programas nunca foram muito complicados.

Os estudantes podem usar esses programas como modelos para escrever os próprios programas no Minitab. Na verdade, a linguagem é tão simples que achamos que qualquer pessoa que programe em outra linguagem consegue ler esses programas e usá-los como modelos do mesmo jeito. Basta pensar nos símbolos c1, c2 etc. como vetores (arrays), em que o \(i\)-ésimo elemento do vetor c1 é acessado por c1(i). Além disso, há as constantes k1, k2 etc.

Um programa do Minitab se chama macro e precisa começar com a instrução gmacro e terminar com a instrução endmacro. A primeira instrução depois de gmacro dá um nome ao programa. Comentários num programa, postos ali para explicar o código, começam com note.

Se o arquivo que contém o programa se chama prog.txt e está gravado no diretório raiz de uma unidade de disco chamada c, então o comando do Minitab

executa o programa. Qualquer saída será impressa na janela Session (se você tiver usado um comando print) ou guardada na planilha (worksheet) do Minitab.

Mais detalhes sobre o Minitab podem ser encontrados no Help do programa. A seguir, damos alguns exemplos de macros do Minitab usadas no texto.

gmacro
bootstrapping
base 34256734
note - original sample is stored in c1
note - bootstrap sample is placed in c2 with each one
note   overwritten
note - medians of bootstrap samples are stored in c3
note - k1 = size of data set (and bootstrap samples)
let k1=15
do k2=1:1000
sample 15 c1 c2;
replace.
let c3(k2)=median(c2)
enddo
note - k3 equals (6.4.5)
let k3=(stdev(c3))**2
print k3
endmacro

gmacro
normalpost
note - the base command sets the seed for the random
note   numbers
base 34256734
note - the parameters of the posterior
note - k1 = first parameter of the gamma distribution
note      = (alpha_0 + n/2)
let k1=9.5
note - k2 = 1/beta
let k2=1/77.578
note - k3 = posterior mean of mu
let k3=5.161
note - k4 = (n + 1/(tau_0 squared) )^(-1)
let k4=1/15.5
note - main loop
note - c3 contains generated value of sigma**2
note - c4 contains generated value of mu
note - c5 contains generated value of coefficient of
         variation
do k5=1:10000
random 1 c1;
gamma k1 k2.
let c3(k5)=1/c1(1)
let k6=sqrt(k4/c1(1))
random 1 c2;
normal k3 k6.
let c4(k5)=c2(1)
let c5(k5)=sqrt(c3(k5))/c4(k5)
enddo
endmacro

let c6=c5 le .5
let k1=mean(c6)
let k2=sqrt(k1*(1-k1))/sqrt(10000)
let k3=k1-3*k2
let k4=k1+3*k2
print k1 k2 k3 k4

gmacro
gibbs
base 34256734
note - data sample is stored in c1
note - starting value for mu.
let k1=mean(c1)
note - starting value for sigma**2
let k2=stdev(c1)
let k2=k2**2
note - lambda
let k3=3
note - sample size
let k4=15
note - n/2 + alpha_0 + 1/2
let k5=k4/2 +2+.5
note - mu_0
let k6=4
note - tau_0**2
let k7=2
note - beta_0
let k8=1
let k9=(k3/2+.5)
note - main loop
do k100=1:10000
note - generate the nu_i in c10
do k111=1:15
let k10=.5*(((c1(k111)-k1)**2)/(k2*k3) +1)
let k10=1/k10
random 1 c2;
gamma k9 k10.
let c10(k111)=c2(1)
enddo
note - generate sigma**2 in c20
let c11=c10*((c1-k1)**2)
let k11=.5*sum(c11)/k3+.5*((k1-k6)**2)/k7 +k8
let k11=1/k11
random 1 c2;
gamma k5 k11.
let c20(k100)=1/c2(1)
let k2=1/c2(1)
note - generate mu in c21
let k13=1/(sum(c10)/k3 +1/k7)
let c11=c1*c10/k3
let k14=sum(c11)+k6/k7
let k14=k13*k14
let k13=sqrt(k13*k2)
random 1 c2;
normal k14 k13.
let c21(k100)=c2(1)
let k1=c2(1)
enddo
endmacro

gmacro
batching
note - k2= batch size
let k2=40
note - k4 holds the batch sums
note - c1 contains the data to be batched (10000 data values)
note - c2 will contain the batch means (250 batch means)
do k10=1:10000/40
let k4=0
do k20=0:39
let k3=c1(k10+k20)
let k4=k4+k3
enddo
let k11=floor(k10/k2) +1
let c2(k11)=k4/k2
enddo
endmacro

gmacro
goodnessoffit
base 34256734
note - generated sample is stored in c1
note - residuals are placed in c2
note - value of D(r) are placed in c3
note - k1 = size of data set
let k1=5
do k2=1:10000
random k1 c1
let k3=mean(c1)
let k4=sqrt(k1-1)*stdev(c1)
let c2=((c1-k3)/k4)**2
let c2=loge(c2)
let k5=-sum(c2)/k1
let c3(k2)=k5
enddo
endmacro

gmacro
dirichlet
note - the base command sets the seed for the random
note   number generator (so you can repeat a simulation).
base 34256734
note - here we provide the algorithm for generating from
note   a Dirichlet(k1,k2,k3,k4) distribution.
note - assign the values of the parameters.
let k1=2
let k2=3
let k3=1
let k4=1.5
let k5=K2+k3+k4
let k6=k3+k4
note - generate the sample with i-th sample in i-th row
note   of c2, c3, c4, c5, ....
do k10=1:5
random 1 c1;
beta k1 k5.
let c2(k10)=c1(1)
random 1 c1;
beta k2 k6.
let c3(k10)=(1-c2(k10))*c1(1)
random 1 c1;
beta k3 k4.
let c4(k10)=(1-c2(k10)-c3(k10))*c1(1)
let c5(k10)= 1-c2(k10)-c3(k10)-c4(k10)
enddo
endmacro