Análise de Dados com o Software R:
Métodos Estatísticos, Computacionais e Econométricos

Prof. Adriano Azevedo Filho (azevedofilho@usp.br)

Análise de distribuição

sumário geral | anterior | próximo

Conteúdo do Módulo

 1 - Lendo o arquivo de dados para o módulo (esalq2012mod.csv)
 
 2 - Introdução à análise da distribuição de dados

 3 - Caracterização gráfica da distribuição de dados
     hist, barplot, boxplot, stripchart
 
 4 - Representação gráfica em eixos cartesianos 
     plot
   
 5 - Medidas de distribuição I: dispersão
     sd, var
 
 6 - Medidas de distribuição II: assimetria e curtose 
 
 7 - Distribuição Normal de densidade de frequência "teórica" ou densidade de probabilidade
     dnorm, pnorm
 
 8 - Distribuição Bernoulli e Distribuição Binomial de distribuição teórica de frequências
     dbinom
 
 9 - Função cumulativa e função quantil
     cumsum, dbinom, pbinom, qbinom, dnorm, pnorm, qnorm, dt, pt, qt, density
 
10 - Estrutura de controle de repetições "for" em funções no R

11 - Aplicação em Economia: medidas de concentração e de desigualdade  
     package ineq

1 - Lendo o arquivo “esalq2012mod.csv” de site na internet

Leia o arquivo “esalq2012mod.csv” que é o mesmo arquivo “esalq2012.csv”, com as modificações nos nomes das variáveis e categorias realizadas por procedimentos descritos aqui.

rm(list=ls())   ## apaga (quase tudo) antes de uma nova análise
alunos<-read.csv2('http://ihbs.com.br/html/esalq2012mod.csv')
alunos[1:5,]

##   sex cur  ing ida cid rep pub pes alt tim sat rmat   pg  car      imc
## 1   m   e 2010  20 isp   1   0  57 172   t   2    s nsei nsei 19.26717
## 2   m   e 2010  20 isp   4   0  78 175   p   2    s nint priv 25.46939
## 3   f   e 2010  22 isp   2   5  59 170   s   2    s pgou priv 20.41522
## 4   m   e 2009  22 isp   6   0  78 179   o   2    s nint nsei 24.34381
## 5   m   e 2010  22 isp   3   4  80 180   t   5    s nsei priv 24.69136

names(alunos)

##  [1] "sex"  "cur"  "ing"  "ida"  "cid"  "rep"  "pub"  "pes"  "alt"  "tim" 
## [11] "sat"  "rmat" "pg"   "car"  "imc"

2 - Introdução à análise de distribuição dos dados

O módulo anterior focou na caracterização da centralidade e posição de conjuntos de dados. Este módulo também focará em conceitos aplicados predominantemente às variáveis quantitativas, com foco na caracterização da * distribuição dos dados: envolvendo dispersão, assimetria, curtose, através de gráficos e medidas-resumo apropriadas.

Os conceitos apresentados neste e nos próximos tópicos focam na caracterização empírica da distribuição dos dados. Estas noções são usualmente utilizadas como estimadores da caracterização teórica. Essa distinção, que pode parecer complicada ou confusa no momento, será esclarecida melhor em outros módulos.

3 - Caracterização gráfica da distribuição: dispersão, assimetria e curtose

A forma mais óbvia para se entender as distribuição dos dados é através de gráficos. No próximo tópico serão apresentadas algumas medidas-resumo importantes que caracterizam aspectos da distribuição dos dados.

3.1 Visualização da dispersão, assimetria e curtose com histogramas e gráficos de barras

O histograma a seguir, mostra exemplos de 2 variáveis (de outro conjunto de dados) com distribuição próxima de simétricas. A primeira delas (à esquerda), contudo, os dados tem uma dispersão maior que no segundo caso (à direita), ao redor de um mesmo valor central. A noção de dispersão tenta capturar a idéia de variabilidade dos dados ao redor de um ponto central.

Haverá mais ou menos dispersão, dentro de uma mesma escala, tanto maior ou menor for a variabilidade ao redor do ponto central. No próximo tópico serão especificadas medidas-resumo que servem para quantificar numericamente a noção de dispersão.

Nos últimos 2 gráficos, o intervalo entre classes é uma unidade, e nos dois casos o histograma apresenta a “densidade” de frequências. É importante recordar que \[\mbox{densidade de frequência} \times \mbox{intervalo entre classes} = \mbox{frequência relativa da classe}\].

A seguir, são caracterizadas, graficamente, as distribuição do peso e do ano de ingresso dos alunos (variáveis “alunos\(pes" e "alunos\)ing”):

par(mfrow=c(1,2))
hist(alunos$pes,breaks=10,col=gray(0.9), xlab="peso em kg",ylab="obs",main="Peso")
barplot(table(alunos$ing), col="bisque", xlab="ano de ingresso",ylab="obs",main="Ano de ingresso")

par(mfrow=c(1,1))

Os 2 gráficos indicam que as variáveis tem dispersão e são relativamente assimétricas. No caso do peso a assimetria é leve e positiva ou assimetrica à direita. No caso do ano de ingresso a assimetria é forte e negativa (ou a esquerda).

Também é interessante caracterizar a curtose da distribuição, ou seja, quanto a distribuição é “achatada” ou “pontuda”, sempre com relação a uma distribuição que tem uma aparência próxima da distribuição normal.

par(mfrow=c(1,1))

As distribuições platicúrticas (curtose negativa) são mais achatadas comparadas ao formato de sino da distribuição Normal (curtose zero). Por outro lado, as distribuições leptocúrticas (curtose posiva) tendem a ser mais “pontudas” comparadas à Normal.

3.2 Box plot ou Gráfico de Tuckey

Os 5 números de Tukey, correspondentes ao percentis 0%, 25%, 50%, 75% e 100%, fa podem ser utilizados para composição de um gráfico que encontra muitas aplicações na estatística, denominado box plot ou “gráfico de Tuckey”. Esse gráfico mostra, simultaneamente, a centralidade, dispersão e assimetria, de uma forma conveniente para comparações.

A seguir são apresentados 2 exemplos de box plots, utilizados para descrição das variáveis peso e altura do data frame “alunos”

par(mfrow=c(1,2))
boxplot(alunos$pes,range=0,col="bisque",main="Peso")
boxplot(alunos$alt,range=0,col="bisque",main="Altura")

par(mfrow=c(1,1))

Com a opção “range=0”, os limites do “box plot” indicam o mínimo e máximo valores observados. Os limites do retângulo central, definem os percentis 25% e 75%, e a linha central indica a mediana ou percentil 50%. Assim, aproximadamente, 50% das observações estão localizadas na região especificada pelo retângulo central. Uma ilustração interessante da interpretação do box plot é apresentada aqui.

Se a opção “range=0” não for utilizada, os limites do box plot serão definidos por no máximo, 1,5 vezes a largura do retângulo central (diferença entre os percentis 75% e 25%), de forma a capturar aproximadamente 99% dos dados. O uso desse procedimento pode evidenciar potenciais “outliers” ou valores muito extremos, que podem ser causados por erros de digitação ou transcrição. A seguir os mesmos gráficos são produzidos sem a opção “range=0”: Os gráficos “box plot” são muito convenientes para comparações, análises de efeitos de variáveis e condicionamentos diversos, algo que faremos no próximo módulo.

3.3 Grafico tipo “Stripchart”

É uma modalidade de gráfico em que os valores são apresentados em uma núvem que se espalha sobre os valores encontrados, de forma a indicar claramente a posição de cada observação, sendo também úteis para um exame cuidadoso dos dados e potenciais efeitos de variáveis e condicionamentos.

A seguir um exemplo de “stripcharts”, envolvendo peso e altura no data frame “alunos”

par(mfrow=c(1,2))
stripchart(alunos$alt,method="jitter",pch=16,xlab="altura (cm)")
stripchart(alunos$pes,method="jitter",pch=16,xlab="peso (kg)")

par(mfrow=c(1,1))

A opção “pch=16” indica o tipo de símbolo que deve ser utilizado para na apresentação (no caso uma bola fechada). Os códigos de símbolos, e outras opções usadas em gráficos, são apresentadas aqui. A opção method=“jitter” faz com que os pontos sejam apresentados em uma núvem, em lugar de alinhados exatamente sobre uma linha. Experimente alterar essas opções.

3.4 Gráfico de densidade ou “density plot” e moda de variável contínua

Em algumas situações, queremos fazer histogramas com muitas classes mas o número de observações é pequeno. Ocorrerá, nesses casos, que muitas classes ficarão sem observações e a consequência disso é um histograma “desdentado”, como ocorre na situação abaixo.

par(mfrow=c(1,2))
hist(alunos$pes,ylab="freq",xlab="altura",main="número adequado de classes",xlim=c(40,120),freq=FALSE)
hist(alunos$pes,breaks=30,ylab="freq",xlab="altura",main="muitas classes",xlim=c(40,120),freq=FALSE)

par(mfrow=c(1,1))

Para sumarizar melhor o formato da distribuição, é possível se estimar uma versão “contínua” do histograma, que tenta capturar as densidades de frequências de uma forma geral, sem grande interferência do número de observações, algo que é apresentado no “gráfico de densidade”.

A seguir apresentamos uma forma padrão de obtenção do “gráfico de densidade” de frequências no R. As áreas dentro desse gráfico, entre dois valores do eixo x, indicam uma estimativa da frequência de observações nessa região.

plot(density(alunos$pes),xlab="peso",ylab="freq",main="Gráfico de densidade - peso")

Moda de variável quantitativa contínua a partir do gráfico de densidade

O grafico de densidade possibilita a definição da moda no contexto geral de uma variável quantitativa contínua, como sendo o valor da variável que leva a densidade de frequência ao máximo valor. A seguinte função consegue esse efeito para esse tipo de variável:

modacont<-function(x){
  z<-density(x)
  z$x[which.max(z$y)] ## retorna o argmax
}

Testando a função na variável peso (alunos$pes)

modacont(alunos$pes)

## [1] 64.09545

plot(density(alunos$pes),xlab="peso",ylab="freq",main="Gráfico de densidade - peso")
abline(v=modacont(alunos$pes)) ## coloca uma linha vertical a partir da moda
arrows(72,0.004,67,0.000,length=0.1)
text(75,0.008,"moda")
text(75,0.006,as.character(round(modacont(alunos$pes),1)))

Combinação de gráficos

É possível combinar um histograma e gráfico de frequências, que utiliza a opção de fazer um gráfico dentro da mesma área gráfica (eixos) previamente definida, através do comando “lines”:

hist(alunos$pes,ylab="freq",xlab="peso (kg)",main="Peso",freq=FALSE,col="bisque",xlim=c(20,140))
lines(density(alunos$pes))

3.5 Exercícios recomendados

• Faça um “box plot” da variável “alunos$imc”, sem utilizar a opção “range=0”. Quantas observações aparecem fora dos limites do box plot?

• Qual deve ter sido a razão de ter sido definida a opção “xlim=c(30,140)” na combinação do histograma e gráfico de densidade, especificada na última sequência de comandos?

4 - Representação gráfica em 2 eixos cartesianos (comando plot)

A visualização em eixos cartesianos é algo de frequente interesse em estatística e outras disciplinas. O entendimento do procedimento computacional utilizado para essa finalidade, ilustrado aqui com o R, é muito similar em outros sistemas ou linguagens.

O princípio geral é estabelecer um conjunto de \(n\) pontos \((x_i, y_i)\), \(i=1,\ldots,n\), correspondentes às coordenadas \(x\) e \(y\) e usar o comando plot para colocação dos pontos no gráfico. Esse princípio será ilustrado aqui na caracterização de funções e na definição de gráficos de dispersão.

4.1 Caracterização de funções \(y=f(x)\) em 2 eixos cartesianos

O caso mais usual é o que envolve 2 dimensões, ou seja, funções do tipo: * \(\displaystyle y=f(x)\) dentro de um certo domínio de interesse com relação a \(x\). É possível também a caracterização em 3 dimensões, mas isso será no momento deixado para outro módulo.

O desenho da função parte da caracterização de \(n\) pontos \((x_i, y_i)\), \(i=1,\ldots,n\), definidos convenientemente, associados ao domínio e contra-domínio da função. Esses pontos são colocados em sistema de eixos cartesianos, nas posições apropriadas e ligados através de segmentos de reta.

“Plotando” a função \(y=x^2\)

Vamos ilustrar esse procedimento, desenhando a função \(y=x^2\) no domínio \(x\in[-2, 2]\), usando a função plot do R que é a principal para gráficos dessa natureza.

• Passo 1: defina um conjunto de pontos dentro do domínio de x, convenientemente espaçados

x<-seq(-2,2,0.1)  ## pontos entre -2 e 2 espaçados de 0.1
x

##  [1] -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7
## [15] -0.6 -0.5 -0.4 -0.3 -0.2 -0.1  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7
## [29]  0.8  0.9  1.0  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2.0

• Passo 2: defina o valor de \(y\) para esses pontos, no contra-domínio da função \(y=x^2\):

y<-x^2
y

##  [1] 4.00 3.61 3.24 2.89 2.56 2.25 1.96 1.69 1.44 1.21 1.00 0.81 0.64 0.49
## [15] 0.36 0.25 0.16 0.09 0.04 0.01 0.00 0.01 0.04 0.09 0.16 0.25 0.36 0.49
## [29] 0.64 0.81 1.00 1.21 1.44 1.69 1.96 2.25 2.56 2.89 3.24 3.61 4.00

• Passo 3: “plote” os pontos usando a função “plot”

plot(x,y,xlab="x",ylab="y")

Sem qualquer opção específica, apenas os pontos definidos em \(x\) e \(y\), tomados na ordem, foram “plotados”. Para obter a “ilusão” visual da função, esses pontos precisam ser ligados por segmentos de reta. Isso pode ser feito automaticamente com o uso da opção “type”, no comando plot.

plot(x,y,xlab="x",ylab="y",type="l")
title(expression(paste("Gráfico da função y=",x^2))) ## coloca o título
grid(col=gray(0.7)) ## coloca o grid da cor selecionada

Com um espaçamento mais largo, os segmentos ficariam mais visíveis, e a “mágica” ficaria evidente, como em

x<-seq(-2,2,1) ## pontos espaçados de 1 unidade
y<-x^2
plot(x,y,xlab="x",ylab="y",type="l")
title(expression(paste("Gráfico da função y=",x^2," com poucos pontos"))) ## coloca o título
grid(col=gray(0.7)) ## coloca o grid da cor selecionada

#### 4.2 Gráficos de dispersão envolvendo 2 variáveis quantitativas

Uma outra aplicação da representação gráfica do comportamento de 2 variáveis quantitativas é a chamada de gráfico de dispersão. Nesse caso, os pontos correspondentes às observações \((x_i, y_i)\), \(i=1,\ldots,n\), estão associados às observações que queremos descrever.

Podemos ilustrar esse procedimento, mostrando um gráfico de dispersão da altura e do peso dos indivíduos no data frame “alunos” usando:

plot(alunos$alt,alunos$pes,xlab="altura (cm)",ylab="peso (kg)")
title(expression(paste("Gráfico de dispersão (altura x peso)"))) ## coloca o título
grid(col=gray(0.7)) ## coloca o grid da cor selecionada

Pode ser conveniente, em alguns casos utilizar símbolos, cores e tamanhos apropriados para caracterizar outras dimensões ou variáveis simultaneamente. Esse é um recurso importante para visualização de situações complexas.

O comando “plot” (e outros comandos gráficos do R) comporta opções para modificação de propridades dos símbolos utilizados como: * pch: altera o tipo de símbolo * col: altera a cor do símbolo * cex: altera o tamanho do símbolo Essas opções podem ser definidas em vetores, com tamanho idêntico à das variáveis utilizadas, de forma que cada observação pode ser particularizada. Veja aqui detalhes sobre opções gráficas disponíveis.

No exemplo abaixo estamos mostrando o mesmo diagrama de dispersão anterior, mas com símbolos da cor rosa e cor azul claro indicando mulheres e homens, bolas e quadrados indicando alunos de economia e agronomia, e o tamanho do símbolo sendo dependente do grau de satisfação com o curso:

corsimbolo<-ifelse(alunos$sex=="f","lightpink","lightblue")
tiposimbolo<-ifelse(alunos$cur=="e",16,15) ## 16 é a bola, 15 é o quadrado
tamsimbolo<-2+(alunos$sat-3)
plot(alunos$alt,alunos$pes,xlab="altura (cm)",ylab="peso (kg)",col=corsimbolo,pch=tiposimbolo,cex=tamsimbolo)
title(expression(paste("Gráfico de dispersão (altura x peso)"))) ## coloca o título
grid(col=gray(0.7))

Para se familiarizar com algumas opções, considere o seguinte exercício, para mostrar os símbolos e cores disponíveis (codificadas numericamente).

x<-1:25
y<-1:25
plot(x,y,col=x,pch=x,cex=2) ## altere o valor de cex e veja o resultado
grid(col=gray(0.7))

5 - Medidas de distribuição I: dispersão

A noção de dispersão é caracterizada numericamente através de medidas-resumo indicadas a seguir, que serão exploradas neste tópico. Essas medidas examinadas neste tópico são relacionadas a seguir.

• medidas-resumo de dispersão: amplitude, desvio absoluto médio, intervalo interquartil, variância e desvio padrão, coeficiente de variação

5.1 Amplitude, desvio absoluto médio e intervalo interquartil

Assuma que

• \(x_1, x_2, \ldots, x_n\) são \(n\) observações de uma variável \(x\).

Esse grupo de medidas é caracterizado pelas seguintes definições:

• amplitude: \(\mbox{máximo}(x)-\mbox{mínimo}(x)\)

• desvio absoluto médio: \(\mbox{DAM}\displaystyle =\frac{\sum_{i=1}^n |x_i - c|}{n}\)

• \(c\) é mais comumente a mediana, mas também pode ser a média ou a moda em alguns contextos. Ele indica a média da distância entre cada observação e a medida de centralidade utilizada, que será minimizada quando essa medida for a mediana.

• intervalo interquartil: percentil 75% - percentil 25% da variável \(x\), que corresponde a diferença entre os quartis 3 e 1.

Essas três medidas de dispersão podem ser computadas no R para a variável altura (alunos$alt) por:

max(alunos$alt)-min(alunos$alt)  ## amplitude

## [1] 47

sum(abs(alunos$alt-median(alunos$alt)))/length(alunos$alt) ## desvio absoluto médio

## [1] 6.625

IQR(alunos$alt) ## intervalo interquartil

## [1] 10.25

5.2 Variância, desvio padrão e coeficiente de variação

Assuma que

• \(x_1, x_2, \ldots, x_n\) são \(n\) observações de uma variável \(x\).

Esse grupo de medidas, que envolvem a variância e medidas relacionadas, é caracterizado pelas seguintes definições:

• variância: \(\displaystyle V(x)=\frac{\sum_{i=1}^n (x_i-\bar x)^2}{n-1}\)

• o uso de \(n-1\) no denominador é uma prática tradicional, uma alternativa válida (e até melhor em algumas situações) é utilizar \(n\). Os símbolos \(s^2\) ou \(V(x)\) podem ser utilizados para indicar a variância no contexto descritivo. O símbolo \(\sigma^2\) é usado no contexto teórico. A variância é uma estimativa da média das distâncias ao quadrado entre cada observação e a média. A unidade da variância é definida em unidades ao quadrado do fenômeno original, algo que torna sua interpretação direta um tanto quanto difícil. Usualmente é mais fácil entender a variância por uma medida associada que é o desvio padrão, descrita a seguir.

• desvio padrão: \(\displaystyle dp(x)=\sqrt{\mbox{V}(x)}\). O desvio padrão é medido na mesma unidade do fenômeno. No contexto da descrição de dados pode ser representado pelos símbolos \(s\) ou \(dp(x)\). No contexto teórico é representado por \(\sigma\). É afetado pela unidade usada para medir o fenômeno.

• coeficiente de variação: \(\displaystyle cv(x)=\frac{dp(x)}{\bar x}\)

• o cv é adimensional e tem sempre o mesmo valor, sendo independente da unidade em que o fenômeno é medido. Não será definido quando \(\bar x=0\). Essas três medidas de dispersão podem ser computadas no R para a variável altura (alunos$alt) por:

var(alunos$alt)  ## variância

## [1] 74.90877

sd(alunos$alt)   ## desvio padrão

## [1] 8.654985

sd(alunos$alt)/mean(alunos$alt) ## coeficiente de variação

## [1] 0.05046113

6 - Medidas de distribuição II: assimetria, curtose e momentos centrais

Também caracterizam a distribuição dos dados, juntamente com a dispersão, as medidas-resumo que estimam numericamente a assimetria (skewness em inglês) e curtose (kurtosis em inglês), as quais serão exploradas neste tópico.

Para as definições assuma que

• \(x_1, x_2, \ldots, x_n\) são \(n\) observações de uma variável \(x\).

5.1 Definições formais de possíveis estimadores da assimetria e curtose

• coeficiente de assimetria: \(\displaystyle g_1(x)=\frac{n}{(n-1)(n-2)}\ \frac{\sum_{i=1}^n (x_i - \bar x)^3}{s^3}\) onde \(s\) é o desvio padrão. Valores próximos de zero sugerem simetria. Valores negativos sugerem assimetria a esquerda e valores positivos sugerem assimetria a direita. Observe que:

• \(\displaystyle \lim_{n\rightarrow\infty} g_1(x)=\frac{1}{n}\frac{\sum_{i=1}^n (x_i - \bar x)^3}{s^3}\)

• coeficiente de curtose: \(\displaystyle g_2(x)=\frac{n(n+1)}{(n-1)(n-2)(n-3)}\frac{\sum_{i=1}^n (x_i - \bar x)^4}{s^4}\ \ -\frac{3(n-1)^2}{(n-2)(n-3)}\) onde \(s\) é o desvio padrão. Valores negativos sugerem distribuições platicúrticas. valores positivos sugerem distribuições mais Leptocúrticas. Observe que:

• \(\displaystyle \lim_{n\rightarrow\infty} g_2(x)= \frac{1}{n}\frac{\sum_{i=1}^n (x_i - \bar x)^4}{s^4}\ \ -3\).

É interessante observar, também, que as fórmulas utilizam, na sua definição, o

• momento central (amostral) de ordem k: \(\displaystyle M_k(x)=\frac{\sum_{i=1}^n (x_i-\bar x)^k}{n}\)
• Note que \(M_2\) é uma possível definição da variância e que os coeficientes de assimetria e de curtose dependem, respectivamente, de \(M_3\) e \(M_4\).

6.2 Implementação no R

As funções definidas acima podem ser implementadas pelas seguintes funções no R:

g1<-function(x){  ## coeficiente de assimetria
  n<-length(x)
  s<-sd(x)
  m<-mean(x)
  n/((n-1)*(n-2))*sum((x-m)^3)/s^3
}

g2<-function(x){  ## coeficiente de curtose
  n<-length(x)
  s<-sd(x)
  m<-mean(x)
  (n*(n+1)/((n-1)*(n-2)*(n-3)))*sum((x-m)^4)/s^4-3*(n-1)^2/((n-2)*(n-3))
}

Testando as 2 funções para as variáveis peso e ano de ingresso (alunos\(pes e alunos\)ing), já descritas por gráficos em tópico acima:

cat("g1=",g1(alunos$pes),"g2=",g2(alunos$pes),"\n")

## g1= 1.031436 g2= 1.448888

cat("g1=",g1(alunos$ing),"g2=",g2(alunos$ing),"\n")

## g1= -3.397879 g2= 16.29778

7 - Distribuição Normal de densidade de frequência “teórica” ou densidade de probabilidade

No início do módulo introduzimos a noção de “gráfico de densidade de frequência”. Essa noção é muito importante em estatística e será explorada futuramente, no contexto da amostragem. Neste tópico, a título de ilustração, é apresentada a distribuição Normal de frequência teórica que tem particular interesse em estatística para caracterização de variáveis contínuas, que serão utilizadas futuramente.

A caracterização descrição algébrica e gráfica dessa distribuição é apresentada a seguir:

• distribuição Normal: \(\displaystyle f(x)=\frac{1}{\sqrt{2 \pi} \sigma}\ e^{-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}}\)

• os símbolos \(\mu\) e \(\sigma\) indicam a média e desvio padrão (teóricos), \(\pi\) é o número pi (3,1416..) e “e” é o número de Euler. Veja mais detalhes sobre essa distribuição aqui.
• A função que define a Normal é definida entre \(-\infty\) e \(\infty\).
• É uma distribuição simétrica com coeficiente de curtose igual a zero

• O R tem uma função interna que representa a distribuição Normal, denominada dnorm(x,\(\mu\),\(\sigma\)), onde, na ordem, \(x\) é o argumento e \(\mu\) e \(\sigma\) representam a média e o desvio padrão associados à distribuição.

Poderíamos “plotar” uma Normal com \(\mu=10\) e \(\sigma=2\) usando:

x<-seq(3,17,0.1)
y<-dnorm(x,10,2)
plot(x,y,type="l",xlab="x",ylab="densidade de freq.",main="Distribuição Normal(10,2)")

A figura acima mostra uma caracterização típica de uma distribuição Normal, onde a média coincide com a moda, definindo o centro da distribuição. O desvio-padrão (\(\sigma\)) é sempre a distância entre a média e o ponto de inflexão da distribuição, caracterizado pelo ponto onde a inclinação da curva altera seu comportamento.

Para calcularmos frequências teóricas, também chamadas probabilidades nesse contexto, da variável com distribuição Normal, com \(\mu\) e \(\sigma\) definidos, estar entre 2 limites, por exemplo \(a\) e \(b\), precisaríamos achar a área entre a curva definida pela função que caracteriza a Normal e o eixo x para essa situação. Ou seja, precisaríamos encontrar \[\int_{a}^{b} \frac{1}{\sqrt{2 \pi} \sigma}\ \ e^{-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}}\ dx\] Infelizmente, nesse caso, a integral indefinida da função que caracteriza a distribuição Normal não existe, e essa integral só pode ser calculada utilizando métodos numéricos aproximados (utilizados para calcular a área).

O R implementa uma função que implementa esses métodos, denominada * pnorm(k,\(\mu\),\(\sigma\)) que calcula a frequência acumulada desde \(-\infty\) até o ponto \(k\).

Assim poderíamos achar a frequência acumulada entre os valores 8 e 13, em uma Normal(10,2) usando

pnorm(13,10,2)-pnorm(8,10,2)

## [1] 0.7745375

x<-seq(3,17,0.1)
y<-dnorm(x,10,2)
plot(x,y,type="l",xlab="x",ylab="densidade de freq.",main="Distribuição Normal(10,2)")
lines(c(8,8),c(0,dnorm(8,10,2)))
lines(c(13,13),c(0,dnorm(13,10,2)))
abline(h=0)
arrows(8,0.04,13,0.04,length=0.1,code=3)
text(10.4,0.03,"área")
text(10.4,0.01,as.character(round(pnorm(13,10,2)-pnorm(8,10,2),4)))

Normal padronizada ou Normal estandartizada - Normal(0,1)

O caso particular da Normal em que temos a média teórica \(\mu=0\) e o desvio padrão teórico \(\sigma=1\) é denominado Normal padronizada ou Normal estandartizada. Esse caso é ilustrado a seguir.

7.3 Exercícios recomendados

• Para variáveis com distribuição Normal (teórica) mostre que, aproximadamente,
• 68% das observações estão no intervalo \(\mu \pm\ \sigma\)
  
• 95% das observações estão no intervalo \(\mu \pm\ 2\ \sigma\)
• 99,7% das observações estão no intervalo \(\mu \pm\ 3\ \sigma\)

8 - Distribuição Bernoulli e Distribuição Binomial de distribuição teórica de frequências

Tal como existem distribuições teóricas de densidades de frequências (ou densidades de probabilidade), para variáveis contínuas tal como a distribuição Normal, vista no tópico anterior, também há em estatística o conceito de distribuição (teórica) de frequências relativas, para o caso de variáveis quantitativas discretas.

Neste tópico abordará 2 dessas distribuições. Se \(X\) é uma variável aleatória discreta, essas distribuições terão as seguintes distribuições:

• Distribuição Bernoulli: \(X=1\) com frequência relativa teórica ou probabilidade \(p\) e \(X=0\) com frequência relativa teórica ou probabilidade \(1-p\). Essa distribuição é muito fundamental em estatística como veremos em outros módulos, sendo a base para muitos modelos estatísticos mais complexos.
• Matemáticamente, pode ser expressa por: \(\displaystyle f_X(x)=p^x\ (1-p)^{1-x}\), \(\ x\in\{0,1\}\) e \(p\in[0,1]\)

• Qualquer fenômeno qualitativo com 2 categorias pode ser representado por uma variável Bernoulli. Por exemplo, no resultado do lançamento de 1 moeda, poderíamos caracterizar por \(X=1\) o resultado “cara” e \(X=0\) o resultado “coroa”. Se as chances de “caras” e “coroas” são as mesmas, temos \(p=0,5\). Esse resultado será melhor desenvolvido em outros módulos.

• Distribuição Binomial: \(\displaystyle f_X(x)={n \choose x} p^x (1-p)^{n-x},\ \ \ x\in\{0,1,\ldots,n\}\) onde \(\displaystyle {n \choose x}=\frac{n!}{(n-x)!\ x!}\) e \(p\in[0,1]\)
• a distribuição Binomial depende de 2 “parâmetros” \(n\) e \(p\)
• note que a distribuição Bernoulli é um caso particular da Binomial onde \(n=1\)
• a Binomial representa a soma de \(n\) variáveis Bernoulli, como veremos em outro módulo.

• Há muitos fenômenos que podem ser caracterizados pela distribuição Binomial: por exemplo, suponha que estamos interessados na distribuição de frequências teóricas do número de resultados “cara” em 10 lançamentos de uma moeda. Claramente podemos ter esse número de “caras” variando de 0 a 10. A distribuição de frequência teórica, nesse caso será definido por uma Binomial com \(n=10\) e \(p=0,5\), como veremos em outros módulos.

Para caracterizar a distribuição Binomial no R, podemos utilizar a função dbinom(x,n,p).

8.2 Exemplo: obtenha as frequências teóricas de uma Binomial com n=10 e p=0,5

Essa Binomial(10,5), como veremos em outro módulo, representa, exatamente, a distribuição de frequências teóricas (ou probabilidades) do número de resultados “cara” observados em 10 lançamentos de 1 moeda cujas frequências teóricas ou probabilidades de “cara” ou “coroa” são 1/2 e 1/2, respectivamente.

Podemos obter essas frequências utilizando:

x<-0:10
dbinom(x,10,0.5)

##  [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 0.2050781250
##  [6] 0.2460937500 0.2050781250 0.1171875000 0.0439453125 0.0097656250
## [11] 0.0009765625

O gráfico de frequências, nesse caso, poderia ser obtido por:

plot(x,dbinom(x,10,0.5), type="h",xlab="x",ylab="freq.",main="Binomial(10,1/2)")

as barras no gráfico indicam as frequências teóricas ou probabilidades em cada situação (e não densidades de frequências ou densidades de probabilidades).

8.3 Exercícios recomendados

• Qual seria a probabilidade de se obter exatamente 10 “caras” em 20 lançamentos de 1 uma moeda?

• Faça o gráfico de uma distribuição Binomial com \(n=100\) e \(p=1/2\), com que distribuição já vista o gráfico se parece?

• Qual seria o procedimento no R para ilustrarmos a distribuição de frequências teóricas ou probabilidades de uma variável com distribuição Bernoulli, com \(p=0,25\)?

9 - Função cumulativa e função quantil

Considerando que a variável de interesse é representada por \(X\), considere as seguintes definições:

• Função cumulativa de frequências ou probabilidades (frequências teóricas): é usualmente representada por \(F(x)\) e indica a frequência relativa “acumulada” de observações de uma dada variável situadas no intervalo \((-\infty,x]\). Equivalentemente, seria a frequência de termos \(X\in(-\infty,x]\mbox{ ou, } X\le x\), para qualquer \(x\in (-\infty,\infty)\).

• O nome da variável pode ser colocado como subscrito da função \(F()\) para deixar claro a que variável ela se refere. É importante observar que estamos usando letras maiúsculas (ex. \(X\)) nesse caso para caracterizar o nome variável e letras minusculas (\(x\)) para indicar os valores que essa variável pode assumir.

• Para uma dada variável X, temos \(F_X(x) \in [0,1]\), \(F_X(-\infty)=0\) e \(F_X(\infty)=1\).

• Função quantil: é usualmente representada por \(Q(q)\) ou \(F^{-1}(q)\), e é definida como a função inversa da Função cumulativa. Ou seja, para qualquer \(q\in [0,1]\) a função quantil nos retornará o valor mínimo de \(x\) que leva a \(F(x)\ge q\). Em muitas situações de interesse o quantil \(Q(q)\) é simplesmente o valor \(x^*\) que atende \(F(x^*)=q\).

A função cumulativa \(F(x)\) e a função quantil \(Q(q)\) são noções que valem tanto no contexto empírico (descritivo) quanto no contexto teórico.

Essas duas funções serão discutidas a seguir dentro do contexto de variáveis quantitativas discretas e contínuas, que envolvem algumas diferenças com relação aos procedimentos de obtenção.

9.1 Função cumulativa \(F(x)\) para variáveis quantitativas discretas

Inicialmente, vamos caracterizar a função cumulativa para uma variável quantitativa discreta, tomando para efeito prático a variável alunos\(ing (ano de ingresso) do data frame "alunos". Outro tópico abordará a definição dessa mesma função para uma variável quantitativa contínua (alunos\)alt).

Com vimos na definição, se representarmos por \(X\) a variável alunos$ing (ano de ingresso), para caracterizarmos a função cumulativa \(F_X(x)\), precisamos conhecer para todos os valores de \(x\in (-\infty,\infty)\) as frequências relativas de \(X\le x\).

Para facilitar a construção de \(F_X(x)\), nesse caso, poderíamos partir das frequências relativas para cada valor possível de \(X\), que é uma variável quantitativa discreta, através de

prop.table(table(alunos$ing))

## 
##       2004       2008       2009       2010       2011 
## 0.01785714 0.01785714 0.10714286 0.26785714 0.58928571

• Para \(x=2009,64\), por exemplo teriamos \(F_X(2009,64)\) como sendo a frequência relativa de observações com ano de ingresso igual ou menor que 2009,64 algo que poderíamos encontrar no R através de

sum(alunos$ing<=2009.64)/length(alunos$ing)

## [1] 0.1428571

que é a soma das frequências relativas associadas a 2004, 2005, 2006, 2007, 2008 e 2009.

Para \(x=-10500\), \(x=-12354,5\) e \(x=1900\), não teremos nenhuma observação de \(X\) que atende a \(X\le x\), logo \(F_X(x)=0\) para esses valores de \(x\). De fato, como no conjunto de dados o menor valor para o ano de ingresso é 2004, para todo \(x<2004\) teremos \(F_X(x)=0\). Da mesma forma, temos \(F_X(x)=1\) para \(x>2011\), dado que só há observações de alunos que entraram até 2011.

• A construção de \(F_X(x)\), nesse caso, pode partir das frequências relativas associadas aos valores dos valores possíveis de \(X\), representados por \(x_1 \le x_2 \le \ldots \le x_n\), já ordenados, que representaremos genericamente por \(f_1, f_2,\ldots, f_n\), da forma indicada na tabela a seguir:

\(i\)	\(F_X(x_i)\)
\(1\)	\(f_1\)
\(2\)	\(f_1+f_2=F_X(x_1)+f_2\)
\(\vdots\)	\(\vdots\)
\(n-1\)	\(f_1+f_2+\ldots+f_{n-1}=F_X(x_{n-2})+f_{n-1}\)
\(n\)	\(f_1+f_2+\ldots+f_{n-1}+f_n=F_X(x_{n-1})+f_n\)

• para a variável alunos$ing (ano de ingresso) os valores de \(f_1, f_2,\ldots, f_n\), correspondentes a \(x_1 \le x_2 \le \ldots \le x_n\), já ordenados, serão definidos por:

prop.table(table(alunos$ing))

## 
##       2004       2008       2009       2010       2011 
## 0.01785714 0.01785714 0.10714286 0.26785714 0.58928571

• para produzir os valores de \(F_X(x_i)\) de acordo com as definições da tabela acima, é necessário fazer a acumulação das frequências até \(x_i\), algo que pode ser conseguido facilmente com apoio da função cumsum do R:

cumsum(c(1,2,3,4,5))

## [1]  1  3  6 10 15

cumsum(prop.table(table(alunos$ing)))

##       2004       2008       2009       2010       2011 
## 0.01785714 0.03571429 0.14285714 0.41071429 1.00000000

Com essas informações, a função \(F_X(x)\) pode ser construída dentro do intervalo \((-\infty, \infty)\) correspondente ao seu domínio:

<td>intervalo</td>
    <td>$F_X(x)$</td>
</tr>

    <td>[$-\infty$, 2004)</td>
    <td>0</td>

<td>[2004, 2008)</td>
    <td>0,01786</td>

<td>[2008, 2009)</td>
    <td>0,03571</td>

<tr>
<td>[2009, 2010)</td>
<td>0,14286</td>

<tr>
<td>[2010, 2011)</td>
<td>0,41071</td>

<td>[2011, $\infty$)</td>
<td>1</td>

A representação gráfica dessa função é apresentada a seguir: A figura destaca a descontinuidade que ocorre no final de cada “degrau” da função. O ponto indica o intervalo fechado. Essa é a típica função cumulativa de uma variável quantitativa discreta.

9.2 Função quantil \(Q(q)\) para variáveis quantitativas discretas

A definição geral da função quantil \(Q(p)\) é especificada a seguir:

• Função quantil \(latex Q(p)\) ou \(latex F^{-1} (q)\): função relacionada a uma dada variável que indica para cada \(q\in [0,1]\) o mínimo valor \(x\) que atende a \(F(x)\ge q\). O domínio dessa função coincide com o domínio da variável de interesse. Essa é a definição geral, cuja aplicação é particularmente relevante para o caso de variáveis quantitativas discretas.

Exemplo

No caso apresentado no útimo tópico, relativo à variável alunos$ing (ano de ingresso), por exemplo, o valor de \(Q(0,2)\) associado a essa variável seria 2010. * Por que \(Q(0,2)=2010\)? esse é o mínimo valor de x (no caso a variável alunos$ing) que atende à condição \(F(x)\ge 0,2\). Como não há um valor que atende exatamente a \(F(x)=0,2\), esse é o melhor que pode ser feito para atendimento à definição. Esse valor também poderia ter sido encontrado diretamente através da função quantile:

quantile(alunos$ing,0.2)

##  20% 
## 2010

Essa situação pode ser ilustrada com o gráfico a seguir:

9.3 Função cumulativa e função quantil para variáveis quantitativas discretas no caso teórico

As noções utilizadas para a situação empírica, vista no tópico anterior, se aplicam diretamente para essa situação, observando que nesse caso:

• \(f_1, f_2,\ldots, f_n\) associados aos valores \(x_1, x_2, \ldots, x_n\) da variável \(X\) serão definidos pela distribuição (teórica) de frequências relativas.

Assim, para encontrarmos a função cumulativa de uma Binomial(10,1/2), cujas frequências relativas seriam definidas por

x<-0:10
dbinom(x,10,0.5)

##  [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 0.2050781250
##  [6] 0.2460937500 0.2050781250 0.1171875000 0.0439453125 0.0097656250
## [11] 0.0009765625

bastaria encontrarmos as frequências acumuladas para cada ponto, utilizando

cumsum(dbinom(x,10,0.5))

##  [1] 0.0009765625 0.0107421875 0.0546875000 0.1718750000 0.3769531250
##  [6] 0.6230468750 0.8281250000 0.9453125000 0.9892578125 0.9990234375
## [11] 1.0000000000

e usarmos procedimentos similares aos utilizados no caso anterior, para “plotar” essa função.

Ocorre, contudo, que o R já tem funções específicas tanto para representar a função cumulativa quanto para a função quantil, para essa situação envolvendo a Binomial: * Frequencias relativas (teóricas): dbinom(x,n,p) * Função cumulativa de frequências teóricas ou probabilidades: pbinom(x,n,p) * Função quantil: qbinom(q,n,p)

Como exemplo de aplicação, obtenha a probabilidade de obter 45 caras ou menos em 100 lançamentos de uma moeda que tem probabilidades 1/2 e 1/2 para caras e coroas, respectivamente.

pbinom(45,100,0.5)

## [1] 0.1841008

Qual seria, nessa última situação, o quantil 0,95 dessa Binomial(100,0.5)?

qbinom(0.95,100,0.5)

## [1] 58

9.4 Função cumulativa e função quantil para variáveis quantitativas contínuas no caso teórico

Os conceitos apresentados neste tópico serão focados na distribuição Normal, mas servem para outras distribuições teóricas contínuas (ex. t-Student, Qui-Quadrado, F). Apenas as funções do R para outras funções serão diferentes, as noções gerais são idênticas.

Considere a definição da distribuição Normal, que corresponde a uma distribuição teórica associada a uma variável contínua \(X\), definida anteriormente por:

• distribuição Normal: \(\displaystyle f(x)=\frac{1}{\sqrt{2 \pi} \sigma}\ e^{-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}}\)

Se \(X\) representa a variável com distribuição Normal(\(\mu\),\(\sigma\)), ou seja, com média teórica \(\mu\) e desvio padrão teórico \(\sigma\), para caracterizarmos a função cumulativa \(F_X(x)\), precisamos conhecer para todos os valores de \(x\in (-\infty,\infty)\) as frequências relativas (teóricas) de \(X\le x\) ou \(X\in (-\infty, x]\). A solução teórica para esse problema é obtida a partir da noção de integral:

• \(\displaystyle F_X(x) = \int_{-\infty}^x \frac{1}{\sqrt{2 \pi} \sigma}\ e^{-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}}\ dx\)

Como a integral indefinida não existe para essa função, a integração é realizada a partir de métodos (numéricos) aproximados para cálculo de área. Como já vimos, no R, essa operação é realizada pela função pnorm, ou seja,

• para \(X\) com distribuição Normal(\(\mu\),\(\sigma\)), \(F_X(x)\) é obtida por pnorm(x,\(\mu\),\(\sigma\)).

Por outro lado, a função quantil Q(q), nesse caso contínuo, que em geral especifica o menor \(x\) que atende a \(F_X(x)\ge q\), em função da continuidade da Normal, será definido pelo valor de \(x\) que atende a \(F_X(x)=q\). No R, a função quantil da Normal é operacionalizada pela função qnorm, também obtida por métodos numéricos, ou seja,

• para \(X\) com distribuição Normal(\(\mu\),\(\sigma\)), \(Q(q)\) é obtida por qnorm(q,\(\mu\),\(\sigma\)).

A função cumulativa e a função quantil são ilustradas a seguir para uma Normal(\(\mu\),\(\sigma\)).

A frequência teórica (ou probabilidade) de valores de \(X\) uma distribuição teórica contínua entre \(a\) e \(b\) é conceitualmente definida por

• frequência teórica de \(X\in [a,b]=F_X(b)-F_X(a)\), onde \(F_X(x)\) é a função cumulativa de \(X\).

Por exemplo, se \(X\) tem distribuição Normal(30,4), podemos achar a frequência teórica (probabilidade) de \(X\in [35, 40]\) e o quantil 0,95 dessa distribuição por

pnorm(40,30,4)-pnorm(35,30,4)  ## probabilidade de X estar entre 35 e 40

## [1] 0.09944011

qnorm(0.95,30,4)  ## quantil 0,95

## [1] 36.57941

• No R os nomes das funções que caracterizam a função cumulativa e função quantil sempre iniciam com as letras p e q, respectivamente.

Exemplo: Aplicação dos conceitos na distribuição teórica t-Student

A distribuição teórica t-Student (veja detalhes aqui) é uma distribuição simétrica muito parecida com a distribuição Normal padronizada. O formato da distribuição é definido por um parâmetro que representaremos por gl. Esse parâmetro gl é usualmente é chamado de “graus de liberdade” da distribuição t-Student normalmente está associado ao tamanho de uma amostra utilizada em testes e definições associadas a intervalos de confiança, em outro módulo.

Para uma distribuição t-Student com \(\nu\) graus de liberdade, teremos no R as seguintes funções: * Função de densidade de frequências teóricas: \(f_X(x)\)=dt(x,gl) * Função cumulativa: \(F_X(x)\)=pt(x,gl) * Função quantil: \(Q(q)\)=qt(q,gl)

Tente reproduzir a figura a seguir com 2 distribuições t-Student com 3 e com 10 graus de liberdade (em cinza e vermelho) e compare-as com uma distribuição Normal(0,1) desenhada em preto. O que acontece na medida que o número de graus de liberdade aumenta? (experimente variar o número de graus de liberdade) * Obtenha a frequência teórica de observações de uma variável que tenha distribuição t-Student com 50 graus de liberdade, entre -1 e 1. Obtenha também os quantis 0,95 e 0,99 dessa distribuição, comparando com os mesmos quantis da distribuição Normal(0,1):

pt(1,50)-pt(-1,50)

## [1] 0.6778744

cat("t-Student(50) Q(0,95) = ", qt(0.95,50)," Q(0,99)=",qt(0.99,50),"\n")

## t-Student(50) Q(0,95) =  1.675905  Q(0,99)= 2.403272

cat("Normal(0,1)   Q(0,95) = ", qnorm(0.95,0,1)," Q(0,99)=",qnorm(0.99,0,1),"\n")    

## Normal(0,1)   Q(0,95) =  1.644854  Q(0,99)= 2.326348

9.5 Função cumulativa e função quantil para variáveis quantitativas contínuas (caso empírico)

Do ponto de vista conceitual, a definição da função cumulativa para o caso de variáveis quantitativas contínuas é o mesmo visto para o caso de variáveis quantitativas discretas.

A operacionalização do conceito, contudo, pode seguir dois caminhos. Um deles é similar ao que foi desenvolvido para o contexto de variáveis discretas. O outro, mais técnico, se baseia na estimativa da função de densidade de frequências (contínua) e aplicação de conceitos utilizados para o caso de distribuições teóricas contínuas.

Há, contudo, algumas tecnicalidades na operacionalização do conceito.

Se \(X\) representa a variável contínua alunos$alt (altura), para caracterizarmos a função cumulativa \(F_X(x)\), precisamos conhecer para todos os valores de \(x\in (-\infty,\infty)\) as frequências relativas de \(X\le x\) ou \(X\in (-\infty, x]\).

Para \(x=170\), por exemplo, poderíamos encontrar \(F_X(170)\) por

sum(alunos$alt<=170)/length(alunos$alt)

## [1] 0.4464286

Para caracterizar a função cumulativa, partiremos da lista ordenada dos \(n\) valores observados de \(X\), \(x_1 \le x_2 \le \ldots \le x_n\), encontrando, inicialmente, \(f_1, f_2, \ldots, f_n\) as frequências relativas associadas a cada observação, usando:

x<-sort(alunos$alt)
f<-x/sum(x)  ## frequencias relativas
f

##  [1] 0.01603332 0.01613743 0.01665799 0.01665799 0.01665799 0.01665799
##  [7] 0.01686622 0.01686622 0.01686622 0.01697033 0.01697033 0.01707444
## [13] 0.01717855 0.01717855 0.01717855 0.01717855 0.01717855 0.01717855
## [19] 0.01738678 0.01759500 0.01769912 0.01769912 0.01769912 0.01769912
## [25] 0.01769912 0.01790734 0.01790734 0.01790734 0.01801145 0.01801145
## [31] 0.01801145 0.01801145 0.01811556 0.01811556 0.01811556 0.01821968
## [37] 0.01821968 0.01821968 0.01821968 0.01821968 0.01821968 0.01821968
## [43] 0.01832379 0.01842790 0.01842790 0.01853201 0.01863613 0.01874024
## [49] 0.01874024 0.01874024 0.01874024 0.01884435 0.01926080 0.01926080
## [55] 0.01978136 0.02092660

y<-cumsum(f) ## frequências cumulativas
y

##  [1] 0.01603332 0.03217074 0.04882874 0.06548673 0.08214472 0.09880271
##  [7] 0.11566892 0.13253514 0.14940135 0.16637168 0.18334201 0.20041645
## [13] 0.21759500 0.23477356 0.25195211 0.26913066 0.28630921 0.30348777
## [19] 0.32087454 0.33846955 0.35616866 0.37386778 0.39156689 0.40926601
## [25] 0.42696512 0.44487246 0.46277980 0.48068714 0.49869859 0.51671005
## [31] 0.53472150 0.55273295 0.57084852 0.58896408 0.60707965 0.62529932
## [37] 0.64351900 0.66173868 0.67995836 0.69817803 0.71639771 0.73461739
## [43] 0.75294118 0.77136908 0.78979698 0.80832900 0.82696512 0.84570536
## [49] 0.86444560 0.88318584 0.90192608 0.92077043 0.94003123 0.95929204
## [55] 0.97907340 1.00000000

A partir dessas definições, poderíamos desenhar um esboço da função cumulativa, usando:

plot(x,y,type="l",xlab="x",ylab="freq. cumulativa",main="Função cumulativa para altura")

Da mesma forma que obtivemos um histograma “desdentado” quando utilizamos muitas classes, aqui temos uma função relativamente “rugosa”, algo que talvez seja inconveniente para efeito de estimativa.

Como procedimento mais técnico, poderíamos utilizar a função de densidade de frequências estimada pela função “density” do R para “melhorar” a estimativa, através dos seguintes procedimentos:

z<-density(alunos$alt)
x<-z$x
f<-z$y
y<-cumsum(f)/sum(f)
lines(x,y,col="red")

A linha vermelha no último gráfico ilustra a nova estimativa, comparada com a situação anterior, onde a “rugosidade” é “alisada”.

10 - Estrutura de controle de repetições “for” em funções no R

Juntamente com o “if” a estrutura de controle “for” é também muito utilizada em programação em geral. O R não é uma exceção.

10.1 Sintaxe do uso do “for”

A estrutura “for” controla a execução de repetições de um certo comando. No R, a sintaxe desse comando é a seguinte:

for (variável in início:fim){
comando1
comando2

comandoN
}

10.2 Exemplos de aplicação do “for”

Execute o exemplo a seguir para ter uma idéia mais concreta de como essa estrutura funciona:

for(i in 1:5) {
    cat("valor de i=",i,"\n")
  }

## valor de i= 1 
## valor de i= 2 
## valor de i= 3 
## valor de i= 4 
## valor de i= 5

Agora um exemplo de uso dentro de uma função. Suponha que deseja criar uma função que reproduza o comportamento do produto vetorial de 2 vetores. Sabemos que dentro do R podemos obter o resultado desejado, usando o operador %*% como a seguir:

x<-c(1,2,3)
y<-c(4,5,6)
x%*%y

##      [,1]
## [1,]   32

Poderíamos contudo, obter o mesmo resultado, através da seguinte função, que assume que os argumentos são vetores de mesmo tamanho:

prodvetor<-function(x,y){
  soma<-0
  n<-length(x)
  for(i in 1:n){
    soma<-soma+x[i]*y[i]
  }
  soma
}

cujo teste, nas mesmas condições anteriores, levaria a

prodvetor(x,y)

## [1] 32

É interessante destacar que em linguagens com o C e o Java, versões similares dessa última função, utilizando o “for”, seriam possibilidades usuais para obtenção do produto vetorial.

11 - Aplicação em Economia: medidas de concentração e de desigualdade

Dentro da Economia (e também em outras disciplinas), há interesse por caracterização de concentração e da desigualdade, em diversas situações. Essas noções, nada mais são que medidas específicas de dispersão, que capturam certos aspectos de interesse da distribuição que, em alguns casos, têm propriedades desejáveis. Algumas das principais medidas serão introduzidas neste tópico, junto com referências para maior aprofundamento sobre o assunto.

• Medidas de concentração: utilizadas para caracterização da intensidade de competição num determinado mercado, a partir do tamanho relativo das empresas que compõe o mercado, ou market-share, medidos convenientemente. Dentro dessas medidas incluem-se: razão de concentração, índice de Herfindahl–Hirschman e outros. Outras disciplinas, como ecologia, por exemplo, usam índices similares para caracterizar a concentração de determinadas espécies, medidos através do índice de diversidade de Simpson.

• Medidas de desigualdade: utilizadas para caracterizar a desigualdade existente na distribuição de um determinado atributo econômico (ex. renda, conhecimento, riqueza etc.) pelos agentes econômicos de um dado contexto. Noções importantes: curva de Lorenz, índice ou coeficiente de Gini, indice de Atkinson, índice de Theil, entre outras.

11.1 Medidas de concentração: razão de concentração \(n\) e índice de Herfindahl–Hirschman

Para efeito da caracterização desse conjunto de medidas, considere 2 mercados em que os vetores \(m1\) e \(m2\) caracterizam a produção das empresas em cada mercado.

m1<-c(55,25,15,5,10,2,1,3,4,5)
m2<-c(12,10,15,10,8,15,12,10,5)

Razão de concentração \(n\) ou RC(n)

A razão de concentração \(n\), representada aqui por RC(n), indica o percentual do mercado representado pelas \(n\) maiores empresas dentro desse mercado. Usualmente utiliza-se 4 ou 8 como valor de \(n\). Veja maiores detalhes aqui.

Se \(P_i\) representam os valores ordenados (em ordem descrescente) de uma medida apropriada da produção da empresa, onde \(i\) e a ordem de classificação, e \(N\) representa o total de empresas desse mercado, temos que:

• \(\displaystyle \mbox{RC}(n) = \frac{\sum_{i=1}^n P_i}{\sum_{i=1}^N P_i}\)

Poderíamos facilmente construir uma função do R para cálculo, usando por exemplo a seguinte função:

RC<-function(x,n){
  sum(sort(x,decreasing=TRUE)[1:n])/sum(x)
}

Ao testar essa função nos mercados \(m1\) e \(m2\) definidos no início deste tópico, obteremos o RC(4) em cada caso por:

RC(m1,4)

## [1] 0.84

RC(m2,4)

## [1] 0.556701

algo que mostra que o mercado \(m1\) é mais concentrado que o \(m2\) pois as quatro maiores empresas representam cerca de 84 % da produção no primeiro caso e cerca de 56% no segundo.

Índice de Herfindahl–Hirschman

O índice de Herfindahl–Hirschman, representado aqui pela sigla IHH, é outra medida, talvez até mais popular que o RC(n), para medida da competitividade em mercados, através de uma medida de concentração.

Se \(f_i\) representa o percentual relativo da produção (ou market-share) de uma empresa \(i\) no mercado como um todo, que tem \(n\) empresas, temos que * \(\displaystyle \mbox{IHH}=\sum_{i=1}^n f_i^2.\)

Se \(f_i\) pode ser medido em taxa ou percentual (ex. 0,10 ou 10%). A interpretação usual (veja mais detalhes aqui) é dada por: * IHH menor ou igual 0,01 (ou 100 na medida percentual): mercado altamente competitivo * IHH entre 0,01 e 0,15 (ou entre 100 e 1500 na medida percentual): mercado desconcentrado * IHH entre 0,15 e 0,25 (ou entre 1500 e 2500 na medida percentual): mercado moderadamente concentrado * IHH acima de 0,25 (ou acima de 2500): mercado altamente concentrado

Implementação no R de função para cálculo do IHH:

IHH<-function(x){
  sum((x/sum(x))^2)
}

Ao testar essa função nos mercados \(m1\) e \(m2\) definidos no início deste tópico, obteremos o IHH em cada caso por:

IHH(m1)

## [1] 0.25952

IHH(m2)

## [1] 0.1197789

Isso mostra que o mercado \(m1\) é altamente concentrado e o mercado \(m2\) é desconcentrado, pela interpretação do IHH apresentada anteriormente. Para obter o IHH na forma “percentual” basta multiplicar esses valores por 10000.

11.2 Medidas de desigualdade: curva de Lorenz, índice de Gini e outros

Há várias formas de caracterizar a desigualdade em economia. Para tanto, partiremos de um conjunto de dados hipotético sobre a renda de 20 indivíduos.

y<-c(9,3,1,3.5,1,0.8,13,1,1,1.5,2,2,3,4,5,6,0.5,3,20,1)

Tal como em situções observadas frequentemente, muitos desses indivíduos tem baixa renda e poucos tem renda muito alta. Ainda que não seja usual, medidas de dispersão já descritas anteriormente poderiam ser utilizadas para caracterizar essa situação. Uma possível descrição gráfica do fenômeno poderia utilizar um histograma como:

hist(y,freq=FALSE,xlab="renda",ylab="freq",main="Renda")

Contudo, é mais frequente a caracterização desse fenômeno através de noções como a curva de Lorenz e índice de Gini, descritos no próximo tópico.

Linha da igualdade de renda, Curva de Lorenz e índice de Gini

Suponha \(r_1\), \(r_2\),\(\ldots\),\(r_n\), representa renda de \(n\) indivíduos, já ordenada em ordem crescente, de forma que

• \(r_1\le r_2 \le \ldots \le r_n\) e

• \(R=\sum_{i=1}^n r_i\) é a renda total.

Com essa notação definimos \(f_k\) e \(x_k\), respectivamente,

• frequência relativa cumulativa da renda até a renda do indivíduo k: \(\displaystyle f_k = \frac{\sum_{i=1}^k r_i}{R}\) para \(k=1,\dots,n\), e \(f_k=0\) para \(k=0\);

• frequência relativa cumulativa de indivíduos até o indivíduo k:
  \(\displaystyle x_k = \frac{k}{n}\), para \(k=1,\dots,n\) e \(x_k=0\) para \(k=0\);

• curva de Lorenz: curva aproximada pela ligação dos pontos \((x_k, f_k)\), para \(k=0,\ldots, n\), como se tivessemos “infinitos” pontos (ou a própria distribuição de probabilidade da renda, numa derivação teórica),

• linha de igualdade de renda: é a curva de Lorenz da situação em que a renda total \(R\) é distribuída equitativamente entre os \(n\) indivíduos, de forma que teremos \(f_k=x_k\), para \(k=0,\ldots,n\).

Para facilitar o entendimento desse conceitos, que são complexos, caracterizaremos inicialmente, a linha de igualdade de renda. Se a renda total \(R\) fosse distribuída igualmente entre os indivíduos, teríamos \(f_k=x_k\), para \(k=0,\ldots,n\). Se representarmos os pontos \((x_k, f_k)\) num gráfico, esses pontos ficarão alinhados sobre uma reta, que passa sobre os pontos (0,0) e (1,1), como na figura abaixo. Se interligarmos esses pontos, chegaremos a figura abaixo, que mostra, graficamente, a chamada “linha da igualdade de renda”, que é a curva de Lorenz da situação mais extrema com relação à igualdade de renda. Nesse caso, os z% mais pobres terão exatamente z% da renda total, para z variando de 0% a 100%. No início do tópico definimos em \(y\) as rendas de 20 indivíduos. Poderiamos, nesse caso, encontrar as frequências relativas cumulativas das rendas ordenadas, através dos seguintes passos:

r<-sort(y) ## ordenando em ordem crescente
f<-cumsum(r)/sum(r) ## achando as frequências relativas cumulativas (renda)
f

##  [1] 0.006150062 0.015990160 0.028290283 0.040590406 0.052890529
##  [6] 0.065190652 0.077490775 0.095940959 0.120541205 0.145141451
## [11] 0.182041820 0.218942189 0.255842558 0.298892989 0.348093481
## [16] 0.409594096 0.483394834 0.594095941 0.753997540 1.000000000

Os dados mostram que a renda dos 4 indivíduos mais “pobres” representa pouco mais de 4% da renda total, e que a renda do 10 indivíduos mais pobres representa cerca de 14,5% da renda total.

Podemos obter a frequência cumulativa dos indivíduos, através de

x<-rep(1,length(f))
x<-cumsum(x)/sum(x) 
x 

##  [1] 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70
## [15] 0.75 0.80 0.85 0.90 0.95 1.00

Isso indica que os primeiros 4 indivíduos representam 20% dos indivíduos e os primeiros 10 indivíduos representam 50% de todos os indivíduos.

Temos agora todos os elementos para caracterizar a curva de Lorenz. Antes disso vamos introduzir um zero no início de cada vetor x e f, para atender à convenção, mencionada na definição, de termos \(x_k=0\) e \(f_k=0\), quando \(k=0\):

x<-c(0,x) ## frequência cumulativa de indivíduos iniciando no zero
f<-c(0,f) ## frequência cumulativa de renda iniciando no zero

A partir dessas definições, colocaremos em um gráfico os pontos \((x_k, f_k)\), representados nos vetores \(x\) e \(f\). Nas figuras a seguir, eixo \(x\) indica a frequência relativa cumulativa de indivíduos e o eixo \(y\) indica a frequência cumulativa da renda (\(f\)). Na figura à esquerda, somente os pontos são representados, na figura à direita, os mesmos pontos são interconectados com segmentos de reta, de forma a aproximar/estimar a chamada curva de Lorenz, localizada abaixo da linha de igualdade de renda. Essas definições permitem a a caracterização do

• Índice de Gini: representado aqui por \(G\), é definido, conceitualmente por \(G = \frac{A}{A+B}\), onde \(A\) é a área delimitada entre a linha da igualdade de renda e a curva de Lorenz e \(B\) é a área abaixo da curva de Lorenz, delimitada pelos eixos, indicadas na figura (à direita) recém apresentada. Como, pela definição da figura, \(A+B=0,5\), temos, como definições alternativas e equivalentes para o índice de Gini: \(G = 2 A\) e \(G = 1-2 B\).

• Limites para o índice de Gini: o índice varia entre 0 e 1, onde 0 indicaria a situação de máxima igualdade de distribuição, com a curva de Lorenz coincidente com a linha de igualdade de renda, e 1 indicaria a situação de máxima desigualdade, com a curva de Lorenz coincidindo com o eixos (n-1 indivíduos com renda zero e 1 indivíduo com toda a renda).

• Valores do índice de Gini variam muito entre países (veja aqui): entre 0,25 na Dinamarca e 0,71 na Namibia, com o Brasil apresentando valor 0,52 em 2012, contra 0,61 em 1999, nessa mesma fonte. Obs: o índice pode ser apresentado na base 0 a 100 (0,1 é reportado como 10), como na fonte citada.

Nota: os comandos utilizados para elaboração dos 2 gráficos apresentados acima são descritos a seguir. Esses comandos utilizam várias técnicas importantes para apresentação de figuras mais elaboradas.

par(mfrow=c(1,2))
plot(x,f,xlab=expression(paste("freq. cumul. de indivíduos - ",x[k])),ylab=expression(paste("freq. cumul. de renda - ",f[k])))
text(0.4,0.7,expression(paste("Pontos (", x[k],", ", f[k],")")))
text(0.4,0.6,"para k = 0,...,n")
plot(x,f,xlab=expression(paste("freq. cumul. de indivíduos - ",x[k])),type="l",ylab=expression(paste("freq. cumul. de renda - ",f[k])),main="Curva de Lorenz e índice de Gini",lwd=2)
lines(c(0,1),c(0,1),lwd=2)
arrows(0.62,0.84,0.72,0.768,length=0.1)
text(0.45,0.90,"linha de igualdade de renda")
abline(h=0)
abline(v=1)
arrows(0.72,0.15,0.68,0.25,length=0.1)
text(0.72,0.08,"curva de Lorenz")
text(0.6,0.4,"área A")
text(0.87,0.3,"área B")
text(0.2,0.60,"índice de Gini")
text(0.2,0.5,"A")
text(0.2,0.4,"A+B")
lines(c(0.15,0.25),c(0.45,0.45))

Procedimentos práticos para computação do índice de Gini

Os últimos parágrafos mostraram os conceitos envolvidos na definição do índice de Gini, a partir da curva de Lorenz. O cálculo das áreas pode ser trabalhoso e sujeito a erros de aproximação. Há muitos procedimentos mais práticos voltados ao cômputo do índice de Gini, por fórmulas que se baseiam nos valores de \(x_k\) e \(f_k\) definidos anteriormente, para \(k=0,\ldots,n\), as frequências relativas cumulativas de indivíduos e renda, respectivamente.

Um dos procedimentos usuais para obtenção do índice de Gini é o que utiliza a * Fórmula de Brown: \(G = 1-\sum_{k=1}^n (x_k-x_{k-1})(f_k+f_{k-1})\).

Essa fórmula usa uma aproximação das áreas pelo método trapezoidal, produzindo um leve erro (para mais) no valor estimado do índice de Gini. Esse erro tenderá a zero na medida que o número de observações aumenta.

Uma possível implementação dessa fórmula no R é apresentada a seguir, com o apoio da estrutura de controle “for”. Note que no R os índices dos vetores sempre iniciam em 1, assim para operacionalizarmos a noção de um vetor com índices variando de 0 a \(n\), utilizamos vetores no R com índices variando de 1 a \(n+1\).

Gini<-function(y){
  n<-length(y)
  f<-c(0,cumsum(sort(y)/sum(y)))
  x<-c(0,cumsum(rep(1,length(y))/length(y)))
  G<-0
  for(k in 2:(n+1)){
    G<-G+(x[k]-x[k-1])*(f[k]+f[k-1])
  }
  G<-1-G
  G
}

Testando a função com o vetor de rendas de 20 indivíduos, definido previamente em \(y\), temos

Gini(y)

## [1] 0.5306888

Há situações, contudo, que os dados de renda são apresentados por classe, onde se sabe a frequência da classe, como na tabela abaixo

Frequência % | Faixa de renda
-------------|----------------
  40         | [0, 1000)
  30         | [1000, 5000)
  20         | [5000, 10000)
   5         | [10000,20000)
   5         | [20000,30000]

Essa é uma situação mais complexa que exigirá procedimentos mais elaborados, descritos na literatura especializada. De um modo geral, estimativas da desigualdade envolvendo a curva de Lorenz e o índice de Gini podem envolver muitas tecnicalidades em função de peculiaridades da variável renda. A própria definição da medida adequada da renda para efeito do cálculo do índice de desigualdade pode ser algo complicado.

11.3 Medidas de desigualdade II: índice de Atkinson e índice de Theil

O índice de Gini é uma medida popular de desigualdade, especialmente no contexto da renda. A literatura nessa área também inclui o indice de Atkinson e o índice de Theil, entre outras. O índice de Theil é interessante em função de algumas propriedades assim como sua relação com o conceito de entropia. Não entraremos mais detalhe sobre esses índices, que podem ser examinados nos links oferecidos e nas referências neles citadas. Uma comparação das várias medidas é apresentada aqui.

O cálculo desses e outros índices pode ser facilitada pela package ineq do R, descrito a seguir.

11.4 Package ineq para índices de concentração e desigualdade

O package ineq do R oferece funções que calculam diretamente a maioria dos indicadores importantes de concentração (IHH e Rosemblunt) , desigualdade (Gini, Atkinson, Theil e outros) e pobreza (Watts, Sen, SST e Foster, não discutidos nesse material). Instale e carregue o package e verifique a ajuda dos principais grupos de indicadores, usando:

require(ineq) ## instale o package antes
?ineq  ## ajuda nas medidas de desigualdade
?conc  ## ajuda nas medidas de concentração
?pov   ## ajuda nas medidas de pobreza

A seguir alguns exemplos de aplicação dos recursos do package para obtenção de medidas de concentração e desigualdade já obtidas anteriormente neste módulo:

require(ineq)
conc(m1,type="Herfindahl")

## [1] 0.25952

ineq(y,type="Gini")

## [1] 0.5306888

ineq(y,type="Atkinson")

## [1] 0.2255963

11.5 Exercício recomendado

• Considerando o vetor \(y\) com rendas de 20 indivíduos, verifique o efeito no índice de Gini de:
• redistribuirmos 2 unidades monetárias de cada um dos 3 indivíduos mais ricos para os 3 indivíduos mais pobres
• aumentarmos a renda de todos eles de 1 unidade
• aumentarmos 10% a renda de cada um
• introduzirmos um imposto de 20% nos indivíduos com renda acima de 10 que seja transformado em benefícios equitativos para todos. Faça suas considerações.