GET00130 - Métodos Computacionais para Estatística II

Suponha que possuímos dados referentes a vendas de produtos em lojas espalhadas por diversas cidades brasileiras.

Na base Vendas.csv possuímos dados do número de um conjunto de produtos vendidos nos meses de Janeiro a Junho, do número de funcionários, da cidade na qual a loja está localizada e se a loja participou de uma campanha publicitária para esses produtos. Todas as informações estão disponíveis para os anos de 2019 e 2020. Caso o número de funcionários seja 0, considere essa informação como NA.

Suponha que desejamos criar um arquivo único com todos os dados disponíveis.

Atividade: Importe o arquivo Vendas.csv, fazendo com que a função de importação leia a variável cod_loja como um character. Armazene em um objeto chamado base.

#Visualizando o objeto base
base

# A tibble: 400 x 11
   cod_loja cidade        campanha   ano venda_JAN venda_FEV venda_MAR venda_ABR
   <chr>    <chr>            <dbl> <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
 1 5654     Niterói              0  2019        73        96        48        59
 2 5654     Niterói              1  2020        76       103        56        61
 3 23029    Belo Horizon…        0  2019        78       104        58        61
 4 23029    Belo Horizon…        1  2020        79       104        59        62
 5 26994    São Gonçalo          0  2019        79       107        61        62
 6 26994    São Gonçalo          0  2020        80       107        63        62
 7 23507    São Gonçalo          1  2019        80       109        63        64
 8 23507    São Gonçalo          0  2020        82       110        63        64
 9 38465    Cabo Frio            1  2019        82       111        64        64
10 38465    Cabo Frio            0  2020        83       113        64        65
# … with 390 more rows, and 3 more variables: venda_MAI <dbl>, venda_JUN <dbl>,
#   funcionarios <dbl>

Fazendo o tratamento da base.

#transformando em factor as variáveis qualitativas
base$cidade = factor(base$cidade)
base$campanha = factor(base$campanha,
                       labels = c("Não","Sim"))

#transformando o 0 em NA na variável número de funcionários
base$funcionarios = na_if(base$funcionarios, 0)


#Visualizando o objeto
base

# A tibble: 400 x 11
   cod_loja cidade        campanha   ano venda_JAN venda_FEV venda_MAR venda_ABR
   <chr>    <fct>         <fct>    <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
 1 5654     Niterói       Não       2019        73        96        48        59
 2 5654     Niterói       Sim       2020        76       103        56        61
 3 23029    Belo Horizon… Não       2019        78       104        58        61
 4 23029    Belo Horizon… Sim       2020        79       104        59        62
 5 26994    São Gonçalo   Não       2019        79       107        61        62
 6 26994    São Gonçalo   Não       2020        80       107        63        62
 7 23507    São Gonçalo   Sim       2019        80       109        63        64
 8 23507    São Gonçalo   Não       2020        82       110        63        64
 9 38465    Cabo Frio     Sim       2019        82       111        64        64
10 38465    Cabo Frio     Não       2020        83       113        64        65
# … with 390 more rows, and 3 more variables: venda_MAI <dbl>, venda_JUN <dbl>,
#   funcionarios <dbl>

Suponha que nosso desejo seja fazer um boxplot da venda no mês de Junho.

#boxplot da variável venda_JUN
base |> 
  ggplot(mapping = aes(y = venda_JUN)) + 
  geom_boxplot() +
  labs(y = "Número de produtos vendidos")  

Ao visualizar os resultados, alguém nos disse que seria interessante comparar as vendas no mês de junho para as lojas que participaram da campanha e para as que não participaram.

#boxplot venda em junho por campanha
base |> 
  ggplot(mapping = aes(x = campanha, y = venda_JUN)) + 
  geom_boxplot() +
  labs(x = "Campanha", 
       y = "Número de produtos vendidos")

Suponha que queiramos observar graficamente a relação entre as vendas dos meses de janeiro e junho para cada cidade em cada ano.

#boxplot venda em junho por campanha
graf_1 = base |> 
  ggplot(mapping = aes(x = venda_JAN, 
                       y = venda_JUN)) + 
  geom_point() +
  labs(x = "Número de produtos vendidos em Janeiro",
       y = "Número de produtos vendidos em Junho")

#Visulaizando gráfico
graf_1

#Criando o gráfico para as combinações de cidade e ano
graf_1 +
  facet_grid(ano ~ cidade)

Atividade: Faça uma comparação gráfica, por meio de histogramas, do comportamento da venda do número de produtos no primeiro trimestre de 2019 por município.

Suponha que o nosso interesse seja o de criar funções e plotar as mesmas no R. Para exemplificarmos, vamos considerar o cenário a seguir.

Vamos considerar uma amostra aleatória simples, \(X_1, \ldots, X_n\) de uma distribuição Exponencial(\(\lambda\)). A função de verossimilhança é dada por

\[L(\lambda|x_1,\ldots,x_n) = \prod_{i=1}^{n} f(x_i|\lambda) = \lambda e^{-\lambda x_1} \times \ldots \times \lambda e^{-\lambda x_n} = \lambda^ne^{-\lambda\sum_{i=1}^{n}x_i}.\]

#Criando uma função que calcula a função de verossimilhança de uma Exponencial 
f.vero = function(l,amostra){
  val = l^length(amostra)*exp(-l*sum(amostra))
  return(val)
}

Suponha que foi observada a seguinte amostra 0.117203051, 0.028553559, 0.006832052, 0.180076271, 0.003652716, 0.040062848, 0.049557241 e 0.007092685.

Vamos fazer o plot da função de verossimilhança de 0 a 100?!

#Plotando o gráfico da função de verossimilhança
ggplot(data = tibble(z = c(0,100)), 
       mapping = aes(x = z)) +
  stat_function(fun = f.vero,
                args = list(amostra = c(0.117203051, 0.028553559, 0.006832052, 0.180076271, 0.003652716, 0.040062848, 0.049557241, 0.007092685))) +
  labs(y = "Função de verossimilhança",
       x = expression(lambda))

Atividade: Faça o gráfico da função de verossimilhança de uma distribuição Gama(\(\alpha\),2), supondo que foi observada a amostra: 0.0633, 0.8730, 0.4798, 0.7240 e 1.5393

Muitas vezes, nos deparamos com resultados teóricos, que podem ser difíceis de visualizarmos. Uma abordagem, que pode nos auxiliar nestes momentos, é tentar verificar esses resultados numericamente.

Para simbolizar a ideia global, suponha que definida uma variável aleatória \(X\) e uma transformação sobre a mesma, suponha \(Y = X^2\). Se existe um resultado que garante qual a distribuição de \(Y\), é possível verificar a mesma de forma numérica da seguinte forma:

• 1 - Gere uma amostra \(x_1,\ldots,x_n\) de \(X\);
• 2 - Aplique a transforação \(Y\) na amostra \(x_1,\ldots,x_n\) obtida no passo anterior;
• 3 - Analise o comportamento da amostra \(y_1,\ldots,y_n\), comparando com a distribuição de \(Y\).

Vamo supor que \(X \sim N(0,1)\) e \(Z \sim N(0,1)\). Existe um resultado que garante que \(Y = X^2 + Z^2 \sim \chi^2_{(2)}\). É possível checarmos o resultado usando o conhecimento de probabilidade, mas aqui, vamos fazer isso numericamente.

# 1 - Gerando n valores de N(0,1)
amostra1 = rnorm(n = 500,
                mean = 0,
                sd = 1)

amostra2 = rnorm(n = 500,
                mean = 0,
                sd = 1)

# 2 - Calculando Y na amostra de X
y = amostra1^2 + amostra2^2
  
# 3 - Analisando o comportamento da amostra e comparando com a distribuição Qui-quadrado com 1 grau de liberdade  
ggplot(data = tibble(valor = y), 
       mapping = aes(x = valor)) +
  geom_histogram(aes(y = ..density..),
                 boundary = 0) +
  stat_function(fun = dchisq,
                args = list(df = 2),
                colour = "Red") +
  labs(x = expression(x^2),
       y = "Densidade")

Atividade: Seja \(X \sim Gama(5,1)\) e \(Y \sim Gama(10,1)\). Verifique numericamente se a distribuição de \(X + Y \sim Gama(15,1)\)

Vamos tentar checar outro resultado?!

Seja \(X_1, \ldots, X_n\) uma amostra aleatória de uma N(\(\mu\),\(\sigma^2\)). Existe um resultado que nos diz que \[md(X_1, \ldots, X_n) \sim N\left(\mu,\frac{\pi\sigma^2}{2n}\right).\]

Vamos verificar o resultado numericamente?! Suponha \(\mu = 10\), \(\sigma^2 = 16\) e \(n = 20\).

#Definindo o número de amostras de tamanho n que serão sorteadas
num.amostra = 5000

#Definindo o tamanho da amostra
n = 20

#Definindo o objeto onde serão armazenads as amostras
amostra = matrix(NA,ncol = num.amostra, nrow = n)

  
#Sorteando as 1000 amostras de tamanho 3
for(i in 1:num.amostra){
  amostra[,i] = rnorm(n = n,
                      mean = 10,
                      sd = 4)
}
  

#Calculando a mediana de cada amostra
variancia = apply(X = amostra,
                  MARGIN = 2,
                  FUN = median)

#Vamos plotar o histograma dos valores observados de T e acrescentar sob o histograma
#a densidade teórica de uma distribuição t
ggplot(data = tibble(variancia), 
       mapping = aes(x = variancia)) +
  geom_histogram(aes(y = ..density..)) +
  stat_function(fun = dnorm, 
                args = list(mean = 10,
                            sd = sqrt(pi*16/(2*n))), 
                col = "red") +
  labs(x = "md",
       y = "Densidade")