import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import warnings
warnings.filterwarnings('ignore')
df_salarios = pd.read_csv('dados/dados_tratados.csv').drop('Unnamed: 0', axis=1)
df_b3 = pd.read_csv('dados/b3_stocks_1994_2020.csv')
df_b3['datetime'] = pd.to_datetime(df_b3['datetime'])
Vimos no artigo anterior que a inferência estatística tem por objetivo fazer generalizações sobre uma população, com base nos dados de uma amostra. Vimos que dois problemas básicos nesse processo são:
Lembrando que parâmetros são funções de valores populacionais, enquanto estatísticas são funções de valores amostrais.
Vamos ver um exemplo.
Exemplo 1:
Uma amostra de $n = 1000$ de um município é escolhida, e cada pessoa da amostra é feita uma pergunta sobre um projeto da prefeitura sobre um determinado problema. A resposta à pergunta pode ser Sim (favorável à solução) ou Não (contrária à solução). O objetivo é estimar a proporção de pessoas favoráveis ao projeto.
Se 700 pessoas responderam Sim à pergunta, a estimativa natural é que a proporção é de $70\%$ ($\frac{700}{1000}$). Essa estimação é baseada na suposição que a amostra é representativa da população. Sabemos que outra amostra poderia levar a outra estimativa. Conhecer as propriedades desses estimadores é um dos propósitos mais importantes da inferência estatística. Vamos ver o que pode ser feito nesse caso particular.
Vamos definir as variáveis aleatórias $X_1, X_2,..., X_n$, tais que:
$$X_1 = \begin{cases} 1, \text{se a i-ésima pessoa na amostra responder Sim}\\ 0, \text{se a i-ésima pessoa na amostra responder Não} \end{cases}$$Sendo $p = P(sucesso)$, onde sucesso significa ser favorável à solução planejada pela prefeitura.
Sabendo que temos uma distribuição binomial com parâmetros $n$ e $p$, e o problema consiste em estimar $p$, temos
$$\hat{p} = \frac{\sum_{i = 1}^n X_i}{n}$$Onde o $\sum_{i = 1}^n X_i$ será a soma dos números de sim
(já que o não (zero) não altera a soma), temos:
No (artigo 7)[https://fhfraga.github.io/pages/conceitos/estatistica/7_introducao_inferencia_estatistica.html], vimos que $\hat{p}$ tem distribuição aproximadamente normal, com parâmetros:
$$E(\hat{p}) = p$$$$\text{Var}(\hat{p}) = \frac{p(1-p)}{n}$$Esses resultados nos ajudam a avaliar as qualidades desse estimador. Por exemplo, a esperança indica que o estimador $\hat{p}$, em média, “acerta” $p$. Dizemos que $\hat{p}$ é um estimador não-viesado (ou não-viciado) de $p$. Ou ainda, a variância indica que para amostras grandes, a diferença entre $\hat{p}$ e p tende a ser pequena, pois para $n \to \infty$ , $\text{Var}(\hat{p}) \to 0$. Nesse caso, dizemos que $\hat{p}$ é um estimador consistente de $p$. Observe que essas propriedades são válidas para o estimador no conjunto de todas as amostras que poderiam ser extraídas da população. Para uma particular amostra, $\hat{p}$ pode estar distante de $p$.
Em algumas situações, podemos ter mais de um estimador para um mesmo parâmetro, e desejamos saber qual deles é “melhor”. O julgamento pode ser feito analisando as propriedades desses estimadores. Vamos exemplificar para facilitar o entendimento.
Exemplo 2:
Desejamos comprar um rifle e, após algumas seleções, restaram quatro alternativas. Para saber qual deles é o melhor foi feito um teste com cada rifle, que consistiu em fixá-lo num cavalete, mirar o centro de um alvo e disparar alguns tiros. Vamos ilustrar com uma imagem.
Se fossemos escolher por precisão, teríamos os dois de cima. Se fossemos escolher pela média dos tiros está dentro dos 3 círculos internos, escolheríamos os dois da esquerda. Então, a arma que escolheríamos seria a do canto superior esquerdo, pois reúne os dois critérios selecionados. Mas se escolhêssemos mais um critério, talvez ela não fosse a escolhida. Muitas vezes, a solução deve ser um conjunto de algumas propriedades.
Vemos na figura termos como precisão e acurácia. A precisão mede a proximidade de cada observação da média de todas as observações, então quanto mais juntos as observações tiverem, mais preciso será. A acurácia mede a proximidade de cada observação do valor alvo que se procura atingir, então, nesse caso, quanto mais próximo ao centro do alvo, mais acurado está.
Um estimador $T$ do parâmetro $\theta$ é qualquer função das observações da amostra, ou seja, $T = g(X_1, ..., X_n$).
Segundo essa definição, um estimador é o que chamamos antes de estatística, porém associando-o a um parâmetro populacional.
O problema da estimação é, então, determinar uma função $T = g(X_1, X_2,..., X_n)$ que seja “próxima” de $\theta$, segundo algum critério. O primeiro critério que abordaremos é dado a seguir.
para todo $\theta$
Se esse valor for diferente, diz-se que $T$ é viesado e a diferença $V(T) = E(T) – \theta$ é chamado o viés de $T$.
Vale notar que a esperança de $T$ é calculada sobre a distribuição amostral de $T$, como tratada no no artigo 7.
Estimativa é o valor assumido pelo estimador em uma particular amostra. Assim, no Exemplo 1, $\hat{p}$ é um estimador de $p$, enquanto $60\%$ é uma estimativa de p.
Uma sequência $\{T_n\}$ de estimadores de um parâmetro θ é consistente se, para todo $\epsilon > 0$, $$P\{|T_n – \theta| > \epsilon \} \to 0, n \to \infty,$$
Aqui temos um caso que nos indica que quanto maior for nossa população, mais nosso erro vai diminuir.
Vamos ver um exemplo com Python.
np.random.seed(1)
# media da população
theta = df_salarios['salario_em_dolares']
# media da amostra com n = 50
t1 = df_salarios['salario_em_dolares'].sample(n=50)
#media da amostra com n = 500
t2 = df_salarios['salario_em_dolares'].sample(n=500)
#media da amostra com n = 1000
t3 = df_salarios['salario_em_dolares'].sample(n=1000)
#media da amostra com n = 2000
t4 = df_salarios['salario_em_dolares'].sample(n=2000)
print(f'n = 50: erro = {np.abs(t1.mean() - theta.mean())}')
print(f'n = 500: erro = {np.abs(t2.mean() - theta.mean())}')
print(f'n = 1000: erro = {np.abs(t3.mean() - theta.mean())}')
print(f'n = 2000: erro = {np.abs(t4.mean() - theta.mean())}')
n = 50: erro = 10520.349880159789 n = 500: erro = 4301.604119840224 n = 1000: erro = 2925.0651198402047 n = 2000: erro = 575.8243801597855
Vemos que a média está tendendo a 0 quanto mais aumentamos o tamanho da nossa amostra. Isso significa que a média da nossa amostra está cada vez mais próxima da média da população.
Isso nos diz que se a variância for menor, mais eficiente é a amostra. Vamos ver um exemplo com Python
Vamos considerar que as amostras que vimos no exemplo acima são não-viesadas
print(f'variância t1 = {t1.var()}')
print(f'variância t2 = {t2.var()}')
print(f'variância t3 = {t3.var()}')
print(f'variância t4 = {t4.var()}')
variância t1 = 3819868903.590204 variância t2 = 4033275973.284533 variância t3 = 4128068881.9699454 variância t4 = 3979110572.27985
Como o resultado, podemos dizer que $t_1$ é a amostra mais eficiente.
também pode ser escrito da forma:
$$ \text{EQM}(T; \theta) = \text{Var}(T) +V²,$$Onde:
t1.var() + (t1.mean() - theta.mean())**2
3930546665.1911817
Um dos procedimentos mais usados para obter estimadores é aquele que se baseia no princípio dos mínimos quadrados.
Vamos a um exemplo para entender como funciona.
Exemplo 3:
Um engenheiro está estudando a resistência $Y$ de uma fibra em função de seu diâmetro $X$ e notou que as variáveis são aproximadamente proporcionais, isto é, elas obedecem à relação $$Y \approx \theta X $$ onde $\theta$ é o coeficiente de proporcionalidade. Baseado numa amostra de cinco unidades, ele deseja estimar o parâmetro $\theta$. Ele submeteu tais amostras a mensuração e testes e obter os seguintes resultados
Observação | X | Y |
---|---|---|
1 | 1,2 | 3,9 |
2 | 1,5 | 4,7 |
3 | 1,7 | 5,6 |
4 | 2,0 | 5,7 |
5 | 2,6 | 7,0 |
Com isso ele obteve $\overline{X} = 1,8$ e $\overline{Y} = 5,4$. Jogando na sua relação ele obtém.
$$Y \approx \theta X <=> \theta = \frac{Y}{X} = \frac{5,4}{1,8} = 3$$Agora, Como verificar a qualidade dessa estimativa? Podemos utilizar o modelo $\hat{Y} = 3X$ e ver como esse prevê os valores de $Y$, para os dados valores de $X$, e como são as discrepâncias entre os valores observados e os estimados pelo modelo. Vamos ver isso numa tabela.
X | Y | 3X | Y - 3x | (Y-3X)² |
---|---|---|---|---|
1,2 | 3,9 | 3,6 | 0,3 | 0,09 |
1,5 | 4,7 | 4,5 | 0,2 | 0,04 |
1,7 | 5,6 | 5,1 | 0,5 | 0,25 |
2,0 | 5,7 | 6,0 | -0,2 | 0,04 |
2,6 | 7,0 | 7,8 | -0,8 | 0,64 |
Total | 0 | 1,06 |
$Y - 3X$ é mede a inadequação do modelo para cada observação da amostra. Porém, vemos que o somatório dá $0$, o que não condiz com a verdade, pois temos erros. $(Y - 3X)²$ é uma tentativa de evitar o problema do sinal, fazendo com que todos os valores fiquem positivos. O seu somatório ($\sum_{i = 1}^{5}(Y_i - 3X_i)²$) é uma tentativa de medir o erro quadrático total da amostra. Quanto menor for o erro quadrático total, melhor será a estimativa. Matematicamente, o problema passa a ser o de encontrar o valor de $\theta$ que minimize a função $$S(\theta) = \sum_{i = 1}^{n}(Y_i - \theta X_i)²$$
Derivando essa função em relação a $\theta$ e igualando o resultado a zero, vamos obter uma equação, que quando resolvida o resultado é: $$\hat{\theta}_{MQ} = \frac{\sum_{i = 1}^{n}X_i Y_i}{\sum_{i = 1}^{n}X_i²}$$
Essa solução é chamada de estimador de mínimos quadrados (EMQ) de $\theta$
Voltando ao exemplo, utilizando essa função, vamos ter:
$$\hat{\theta}_{MQ} = \frac{(1,2 \cdot 3,9) + (1,5 \cdot 4,7) + (1,7 \cdot 5,6) + (2,0 \cdot 5,7) + (2,6 \cdot 7,0)}{1,2² + 1,5² + 1,7² + 2,0² + 2,6²}$$$$\hat{\theta}_{MQ} = \frac{50,85}{17,34} = 2,93$$Substituindo pelo valor de $\theta$ correto na tabela
X | Y | 2,93X | Y - 2,93x | (Y-2,93X)² |
---|---|---|---|---|
1,2 | 3,9 | 3,5 | 0,4 | 0,16 |
1,5 | 4,7 | 4,4 | 0,3 | 0,09 |
1,7 | 5,6 | 5,0 | 0,6 | 036 |
2,0 | 5,7 | 5,9 | -0,2 | 0,04 |
2,6 | 7,0 | 7,6 | -0,6 | 0,36 |
Total | 0,5 | 1,01 |
Achando o valor de $S(\theta) = 1,01$, que é menor que o valor anterior ($S(\theta) = 1,06$)
Lógico que neste exemplo, não esperávamos uma relação perfeita entre as duas variáveis, já que o diâmetro da fibra não é o único responsável pela resistência. Porém, estamos supondo que, para um dado valor da variável explicativa $X$, os valores da variável resposta $Y$ seguem uma distribuição de probabilidade $f_Y (y)$, centrada em $\theta X$. Isso equivale a afirmar que, para cada $X$, o desvio $\epsilon = Y – \theta X$ segue uma distribuição centrada no zero. Podemos, então, escrever: $$E(Y|x) = \theta x, \text{para todo valor de } x$$
É comum supor que $\epsilon$ tem a mesma distribuição, para todo valor $x$ da variável explicativa $X$. Desse modo, é comum escrever $$Y = \theta x + \epsilon$$ com $\epsilon$ segundo a distribuição $f_{\epsilon}(.)$, com média zero. Como ilustração, poderíamos supor que $\epsilon \sim N(0, \sigma²)$, para todo $x$. Quanto menor for a variância $\sigma²$, melhor será a “previsão” de $Y$ como função de $x$. Assim, parece razoável escolher $\theta$ que torna mínima a soma dos quadrados dos erros: $$\sum_{i = 1}^n \epsilon_i ^2= \sum_{i = 1}^n (Y_i - \theta X_i)^2$$
O modelo acima pode ser generalizado, de modo a envolver outras funções do parâmetro $\theta$, resultando no modelo $$Y = g(X; \theta) + \epsilon$$ e devemos procurar o valor de $\theta$ que minimize a função $$\sum_{i = 1}^n \epsilon_i ^2= \sum_{i = 1}^n (Y_i - g(X_i; \theta))^2$$
Uma amostra verossímil é uma amostra que fornecesse a melhor informação possível sobre um parâmetro de interesse da população, desconhecido, e que desejamos estimar.
O princípio da verossimilhança afirma que devemos escolher aquele valor do parâmetro desconhecido que maximiza a probabilidade de obter a amostra particular observada, ou seja, o valor que torna aquela amostra a “mais provável”. ”. O uso desse princípio conduz a um método de estimação pelo qual se obtêm os chamados estimadores de máxima verossimilhança que, em geral, têm propriedades muito boas.
Exemplo 4: Vamos supor que temos $n$ pesquisas sobre a intenção de voto de um determinado candidato, onde Sim é sucesso ($P(\text{Sucesso}) = p$, $0 < p < 1$ e $X = $ números de sucessos. Devemos tomar como estimador aquele valor de $p$ que torna a amostra observada mais provável de ocorrer. Vamos supor $n = 3$ e temos dois sucessos e um fracasso. A função de verossimilhança é: $$L(p) = P(\text{2 sucessos e 1 fracasso}) = p^2(1-p) = p² - p³$$
Maximizando essa função em relação a $p$, temos
$$L'(p) = 2p - 3p²$$Achando as raízes da função, temos:
$$x_1 = 0, x_2 = \frac{2}{3}$$Como queremos o ponto de máximo, temos $\hat{p} = \frac{2}{3}$, que é o estimador de verossimilhança (EMV) de p.
Então, podemos definir como
Definição: A função de verossimilhança é definida por $$L(\theta; x_1,...,x_n) = f(x_1;\theta),..., f(x_n;\theta)$$ O estimador de máxima verosimilhança de $\theta$ será o valor de $\theta_{MV}$ que maximiza $L(\theta; x_1,...,x_n)$
Como vimos acima o exemplo da binomial, o EMV do parâmetro $p$ de uma distribuição binomial é $$\hat{p}_{MV} = \frac{X}{n}$$ Onde:
Para uma variável aleatória com distribuição de Bernoulli, o estimador de máxima verossimilhança para $p$ é a proporção amostral (que segue a mesma definição de média amostral): $$\hat{p}_{MV} = \frac{\sum_{i = 1}^n X_i}{n} = \hat{p}$$
Para uma variável com distribuição normal, o EMV para média é $$\hat{mu}_{MV} = \frac{\sum_{i = 1}^n X_i}{n} = \overline{X}$$ Ou seja, para uma população com distribuição normal, a média amostral é o estimador de máxima verossimilhança (EMV) para a média populacional.
E o EMV para a variância, para uma distribuição normal, é: $$\hat{\sigma}²_{MV} = \frac{\sum_{i = 1}^n (X_i - \overline{X})^2}{n}$$
Para uma variável aleatória com distribuição de Poisson, o EMV também é a média amostral: $$\hat{\lambda}_{MV} = \frac{\sum_{i = 1}^n X_i}{n} = \overline{X}$$
Para uma variável com distribuição exponencial, a média é o inverso do parâmetro: $E(X) = \frac{1}{\lambda}$, então a estimativa de máxima verossimilhança para o parâmetro 𝜆 é o inverso da média amostral: $$\hat{\lambda}_{MV} = \frac{1}{\overline{X}}$$
Para uma variável uniforme no intervalo (0, $\theta$), o estimador de máxima verossimilhança para $\theta$ é o maior valor observado na amostra.
Todos os estimadores vistos até agora foram estimadores pontuais. Eles apenas especificam um valor para o estimador. Esse procedimento não permite julgar qual possível magnitude do erro que estamos cometendo. Daí surge os intervalos de confiança, que serão construídos baseados na distribuição amostral do estimador pontual.
Vamos supor que queremos estimar a média $\mu$ de uma população qualquer, e para tanto usarmos a média $\overline{X}$ de uma amosta de tamanho $n$. Do teorema do limite central, $$e = (\overline{X} - \mu) \sim N(0, \sigma²_x),$$ com $\text{Var}(\overline{X}) = \sigma^2 = \frac{\sigma²}{n}$. Daqui podemos determinar qual a probabilidade de cometermos erros de determinadas magnitudes. Por exemplo, $$P(|e| < 1,96\sigma_x) = 0,95$$ ou $$P(|\overline{X} - \mu| < 1,96\sigma_x) = 0,95$$ que é equivalente a $$P(-1,96\sigma_x < \overline{X} - \mu < 1,96\sigma_x) = 0,95$$ e, finalmente, $$P(\overline{X}-1,96\sigma_x < \mu < \overline{X}+1,96\sigma_x) = 0,95$$
Convém lembrar que $\mu$ não é uma variável aleatória e sim um parâmetro, e a expressão acima deve ser interpretada da seguinte maneira: se pudéssemos construir uma quantidade grande de intervalos (aleatórios!) da forma $]\overline{X} – 1,96\sigma⎯x ,\overline{X} + 1,96\sigma_x[$, todos baseados em amostras de tamanho $n$, 95% deles conteriam o parâmetro $\mu$. Dizemos que $\gamma = 0,95$ é o coeficiente de confiança.
Exemplo 5: Uma máquina enche pacotes de café com uma variância igual a $100 g²$. Ela estava regulada para encher os pacotes com 500 g, em média. Agora, ela se desregulou, e queremos saber qual a nova média $\mu$. Uma amostra de 25 pacotes apresentou uma média igual a 485 g. Vamos construir um intervalo de confiança com 95% de confiança para $\mu$
$$\sigma_x = \frac{\sigma}{\sqrt{n}} = \frac{10}{\sqrt{25}} = \frac{10}{5} = 2$$$$IC(\mu; 0,95) = 485 \pm 1,96 \cdot 2$$$$IC(\mu; 0,95) = ]481, 489[$$Então temos que a nova média está entre o intervalos de 481 g e 489 g.
Se $T$ for um estimador do parâmetro $\theta$, e conhecida a distribuição amostral de $T$, sempre será possível achar dois valores $t_1$ e $t_2$, tais que $$P(t_1 < \theta < t_2) = \gamma$$
Para uma dada amostra, teremos dois valores fixos para $t_1$ e $t_2$, e o intervalo de confiança para $\theta$, com coeficiente de confiança $\gamma$, será indicado do seguinte modo: $$IC(\theta, \gamma) = ]t_1, t_2[$$
Se a variância populacional $\sigma²$ não for conhecida, podemos substituir $\sigma_x$ por $\frac{S}{\sqrt{n}}$, então teremos $$P\Big(\overline{X}-1,96\cdot\frac{S}{\sqrt{n}}< \mu < \overline{X}+1,96\cdot\frac{S}{\sqrt{n}}\Big) = 0,95$$ onde S2 é a variância amostral.
Para $n$ grande, da ordem de 100, a equação acima pode ainda ser usada. Para $n$ não muito grande, a distribuição normal não pode mais ser usada e terá de ser substituída pela distribuição t de Student.
Definição: Para um coeficiente de confiança qualquer $\gamma$, teremos de usar o valor $z(\gamma)$ tal que $P(z(\gamma) < Z < z(\gamma)) = \gamma$, com $Z \sim N(0, 1)$. O intervalo fica:
Para média $$IC(\mu, \gamma) = \Big]\overline{X} - z(\gamma)\cdot\frac{\sigma}{\sqrt{n}}; \overline{X} + z(\gamma)\cdot\frac{\sigma}{\sqrt{n}}\Big[$$
Para a proporção temos: $$IC(p, \gamma) = \Big]\hat{p} - z(\gamma)\cdot \sqrt{\frac{p\cdot(1-p)}{n}}; \hat{p} + z(\gamma)\cdot \sqrt{\frac{p\cdot(1-p)}{n}}\Big[$$
Exemplo 6:
Numa pesquisa de mercado, n = 400 pessoas foram entrevistadas sobre determinado produto, e 60% delas preferiram a marca A. Aqui, $\hat{p} = 0,6$ e um intervalo de confiança para $p$ com coeficiente de confiança $\gamma = 0,95$ será:
$$IC(p, \gamma) = \Big]\hat{p} - z(\gamma)\cdot \sqrt{\frac{p\cdot(1-p)}{n}}; \hat{p} + z(\gamma)\cdot \sqrt{\frac{p\cdot(1-p)}{n}}\Big[$$$$IC(p, 0,95) = \Big]0,6 - 1,96 \cdot \sqrt{\frac{0,6\cdot 0,4}{400}}; 0,6 + 1,96 \cdot \sqrt{\frac{0,6\cdot 0,4}{400}}\Big[$$$$IC(p, 0,95) = ]0,55; 0,65[$$Com python teriamos:
def intervalo_de_confianca_para_proporcao(proporcao, coeficiente_confianca, tamanho_amostra):
confianca = coeficiente_confianca * np.sqrt((proporcao * (1-proporcao)) / tamanho_amostra)
return proporcao - confianca, proporcao + confianca
# para nível de confiança de 95%
proporcao = len(df_salarios[df_salarios['nivel_experiencia'] == 'Pleno']) / len(df_salarios)
intervalo_de_confianca_para_proporcao(proporcao, 1.96, len(df_salarios))
(0.2012542903652969, 0.22750736076652733)
Vimos que, obtida a distribuição amostral de um estimador, podíamos calcular a sua variância. Para achar o desvio padrão, basta tirar a raiz quadrada da variância. O desvio padrão de um estimador pode ser chamado de erro padrão. Se $T$ for um estimador do parâmetro $\theta$, chamaremos de erro padrão de $T$ a quantidade
$$EP(T) = \sqrt{\text{Var}(T)} $$O erro padrão (ou desvio padrão) da média amostral é dado por:
$$EP(\overline{X}) = \frac{\sigma}{\sqrt{n}}$$Quanto não conhecemos a variância, podemos oberter o erro padrão estimado da média amostral, que é:
$$\hat{EP}(\overline{X}) = \frac{S}{\sqrt{n}}$$O erro padrão da proporção amostral é dado por:
$$EP(\hat{p}) = \sqrt{\frac{p\cdot(1-p)}{n}}$$Para construir intervalos de confiança e testar hipóteses será necessário conhecer a distribuição amostral dos estimadores. Como só temos um conjunto de dados e não dados hipotéticos, estas distribuições amostrais terão de ser obtidas de outra maneira. Usualmente isso é feito usando teoremas como o Teorema do Limite Central, obtendo-se uma distribuição aproximada para os estimadores, que vale para tamanhos de amostras grandes.
A crítica que se faz à teoria frequentista é a possibilidade de “replicar dados”, bem como o recurso à teoria assintótica. Uma teoria que não faz uso de tais argumentos é a inferência bayesiana.
Na inferência Bayesiana, além das informações disponíveis na amostra, que chamamos de distribuição a posteriori, incorporamos informações que sabemos ser verdadeiras, mas que não estarão representadas na amostra (não observáveis). Essas informações podem ser denominadas subjetivas e compõem a distribuição a priori.
Por exemplo, nas eleições, sabemos que um candidato não terá 100% dos votos. Essa informação pode ser incorporada ao modelo para aprimorar a estimativa. Por outro lado, a inserção de informações que não são verdadeiras torna o modelo inapropriado. Por exemplo, se for considerado que nenhum candidato terá 70% dos votos, mas se isso for de fato possível, a estimativa estará enviesada.
Na inferência Bayesiana, considera-se também que o parâmetro populacional a ser estimado é uma variável aleatória (e não fixo como na inferência clássica). Ademais, o que chamávamos de intervalo de confiança na inferência clássica passa a ser chamado de intervalo de credibilidade.
Vamos ver um exemplo
Exemplo 7: Vamos imaginar que um e-commerce está com um produto novo de uma determinada marca, que chamaremos de produto A. Porém, essa marca não é tão conhecida. Esse e-commerce nos pede para que seja feito um estudo sobre a viabilidade de mandar disparos de CRM sobre esse produto para clientes que compraram um produto similar da marca B. O objetivo é entender se clientes que compraram a marca B se interessariam pela produto da marca A. Vamos ver como faríamos isso com inferência bayesiana.
Primeiro devemos ter os seguintes parâmetros:
Vamos analisar os dados do e-commerce (fictício) e lá vemos o seguinte caso: O produto A tem 987 avaliações, sendo:
O produto B tem 1359 avaliações, sendo:
Vamos considerar 5 e 4 estrelas como avaliação positiva e 2 e 1 estrela como avaliação negativa. 3 estrelas vamos desconsiderar por ser o meio termo, onde não sabemos se é mais positivo ou mais negativo.
# Dados produto A
X_a = np.array([1] * (int(0.84 * 987) + int(0.05 * 987)) + [0] * (int(0.04 * 987) + int(0.02 * 987)))
np.random.shuffle(X_a) # apenas para embaralhar os dados
# Dados produto B
X_b = np.array([1] * (int(0.85 * 1359) + int(0.02 * 1359)) + [0] * (int(0.05 * 1359) + int(0.02 * 1359)))
np.random.shuffle(X_b)
O que o pessoal do e-commerce quer saber é: "Qual a probabilidade das pessoas que compraram o produto B gostarem mais de comprar o produto A". Pra simplificar, imagina que você quer comprar uma TV, você já comprou uma TV nesse e-commerce da marca B. Acabou de lançar uma TV da marca A. Vamos dizer que você já quer trocar de TV. Qual a probabilidade de você gostar mais da TV da marca A do que de uma da marca B?
E isso é fundamental porque a priori é baseado na sua experiência, já que se você já conhece a marca, você tem uma opinião sobre ela, seja positiva ou negativa, enquanto para uma marca que você não conhece, você ainda não tem uma opinião formada.
Se não conhecemos a marca A, temos a mesma probabilidade de avaliar positivamente ou negativamente. Para modelar isso, vamos utilizar a distribuição de probabilidade Beta, que é uma distribuição conveniente para comportamentos aleatórios de porcentagens e proporções. Para escolher os parâmetros, pode usar essa calculadora de distribuição beta. Escolhemos os parâmetros:
Para modelar a posteriori vamos usar a distribuição de probabilidade de Bernoulli, pois a mesma é uma distribuição usada para proporções, onde temos 1 (sucesso) e 0 (fracasso). No nosso caso o 1 será para avaliações positivas e o 0 para avaliações negativas.
with pm.Model() as model:
# priori
priori_a = pm.Beta('produto_a', 1, 1)
priori_b = pm.Beta('produto_b', 5, 2)
# Deterministico
delta = pm.Deterministic('delta', priori_a - priori_b)
# posteriori
obs_a = pm.Bernoulli('obs_a', priori_a, observed=X_a)
obs_b = pm.Bernoulli('obs_b', priori_b, observed=X_b)
# likelihood
trace = pm.sample(draws = 10000, tune=5000, step=pm.Metropolis())
Multiprocess sampling (4 chains in 4 jobs) CompoundStep >Metropolis: [produto_b] >Metropolis: [produto_a]
Sampling 4 chains for 5_000 tune and 10_000 draw iterations (20_000 + 40_000 draws total) took 17 seconds. The number of effective samples is smaller than 25% for some parameters.
pm.plot_posterior(trace['produto_a'])
<Axes: title={'center': 'x'}>
Temos um valor 94% para avaliar o produto A positivamente com um intervalo de credibilidade que vai de 92% a 95%.
pm.plot_posterior(trace['produto_b'])
<Axes: title={'center': 'x'}>
Temos um valor 93% para avaliar o produto B positivamente com um intervalo de credibilidade que vai de 91% a 94%.
plt.figure(figsize=(14, 10))
ax1 = plt.subplot(3, 1, 1)
ax2 = plt.subplot(3, 1, 2)
ax3 = plt.subplot(3, 1, 3)
ax1.hist(trace['produto_a'], histtype='stepfilled', bins=40, density=True)
ax1.set_xlim([0.84, 1])
ax1.vlines(trace['produto_a'].mean(), 0, 60, linestyle='--', color='black')
ax1.set_title('Posteriori produto A')
ax2.hist(trace['produto_b'], histtype='stepfilled', bins=40, density=True)
ax2.set_xlim([0.84, 1])
ax2.vlines(trace['produto_b'].mean(), 0, 60, linestyle='--', color='black')
ax2.set_title('Posteriori produto B')
ax3.hist(trace['delta'], histtype='stepfilled', bins=40, density=True)
ax3.vlines(0, 0, 35, linestyle='--', color='black')
ax3.set_title('Posteriori Delta')
plt.show()
print(f'Probabilidade do cliente gostar mais do produto A: {np.mean(trace["delta"] > 0)}')
print(f'Probabilidade do cliente gostar mais do produto B: {np.mean(trace["delta"] < 0)}')
Probabilidade do cliente gostar mais do produto A: 0.869575 Probabilidade do cliente gostar mais do produto B: 0.130425
Vemos que com a inferência bayesiana, vale a pena mandar disparos de CRM para o cliente que já comprou o produto B, pois ele tem a probabilidade de 86,7% de avaliar positivamente (gostar) do produto A.