Generalized additive model, quando não sabemos que modelo usar.

Vamos olhar algo que um amigo meu me perguntou sempre sobre, a algum tempo atrás, o GAM, ou Generalized additive model.

Primeiro, vamos ver qual o interesse em outro método estatístico. Sabemos que existem muitos modelos matemáticos, que explicam partes determinísticas de processos biológicos, como crescimento populacional com o modelo do Lotka-Volterra ou o comportamento de populações no espaço com o modelo de metapopulações do Levins.

Normalmente a gente não coleta dados e ajusta um modelo de crescimento populacional neles. A gente usa coisas como análise de variância e regressão linear entrou outros modelos, porque são modelos mais simples, mais robustos que acabamos por detectar padrões nos dados, e depois vamos mais afundo nesses padrões.

Vamos pensar num exemplo aqui. Vamos gerar dados de um preditor, e duas respostas diferentes.

1
2
3
4
5
###Simulando Dados
set.seed(123)
preditor <- runif(30,0,4*pi)
resposta_pot<-5+2*(preditor)^2+rnorm(30,0,10)
resposta_seno<-5+2*sin(preditor)+rnorm(30)

Seja o preditor a quantidade de açúcar no néctar e a resposta a territorialidade na fonte desse néctar. Suponha que temos algo assim:

Bem a gente pode ajustar uma reta ali, usando uma regressão linear e vai dar um resultado legal.

E teremos uma relação significativa.

Call: lm(formula = resposta_pot ~ preditor) Residuals: Min 1Q Median 3Q Max -43.001 -24.431 -2.866 23.330 53.420 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -52.67 11.82 -4.456 0.000123 *** preditor 26.04 1.47 17.718 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 28.98 on 28 degrees of freedom Multiple R-squared: 0.9181, Adjusted R-squared: 0.9152 F-statistic: 313.9 on 1 and 28 DF, p-value: < 2.2e-16

Com um valor de R quadrado bastante alto, e fim de papo? Pode ser que sim, mas podemos fazer melhor. Veja que nesse caso, temos um modelo exponencial.

Apesar da regressão resolve ali localmente, o primeiro problema é o fato de que não podemos extrapolar baseado nesse modelo de regressão, já que ele vai falhar miseravelmente além do intervalo de dados que temos, é so pensar que a curva vai subir muito mais “rápido” que a reta da regressão para imaginar isso, outro problema é que nosso residuo não tem uma distribuição normal, temos um padrão forte nele, o padrão que vem do determinismo dos dados que não está sendo coberto.

Na figura Residual vs Fitted a gente vê do que estou falando, aqui era para nao ter padrão nenhum, no entanto temos uma curva bem fácil de identificar até no olhômetro.

Mas esse é quando a gente ta com sorte, agora imagine se os dados tem um comportamento todo estranho, um padrão, mas que não é algo com cara de reta.

Ai a regressão não vai dar em nada, no caso anterior a gente poderia resolver com uma transformação, mas nesse caso não.

E aqui a robustez não funciona

Call: lm(formula = resposta_seno ~ preditor) Residuals: Min 1Q Median 3Q Max -2.8625 -1.3596 -0.2053 1.5895 3.9992 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.46404 0.77687 7.033 1.19e-07 *** preditor -0.06799 0.09658 -0.704 0.487 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.904 on 28 degrees of freedom Multiple R-squared: 0.01739, Adjusted R-squared: -0.0177 F-statistic: 0.4956 on 1 and 28 DF, p-value: 0.4873

Mas o padrão está la.

Denovo a gente vê que não deveria acreditar na regressão, a primeira figura do Residual vs Fitted, a gente vê que estamos quebrando premissas, não há uma distribuição normal nos resíduos, tem um padrão neles, mas não conseguimos pegar no modelo.

Então é pra esse tipo de coisa que existe o GAM, a gente sabe que existe um padrão determinístico nos dados, mas não sabemos qual é ele, então ou a gente sai chutando vários modelos, ou a gente uma o Generalized additive model. De certa forma é bem simples usar, ele tem uma implementação no R no pacote mgcv, que é bem simples de usar.

a sintaxe é praticamente a mesma de um glm ou lm, so adicionando aquela função s(), então o ajuste fica assim.

> modelo_gam<- gam(resposta_pot~s(preditor,k=5,fx=T),family = gaussian) > summary(modelo_gam) Family: gaussian Link function: identity Formula: resposta_pot ~ s(preditor, k = 5, fx = T) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 134.639 1.618 83.22 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(preditor) 4 4 908.1 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.992 Deviance explained = 99.3% GCV = 94.22 Scale est. = 78.516 n = 30

E temos um sumario bem parecido com o da regressão linear, o modelo significativo e tal, e podemos plotar o modelo.

E ele consegue capitar bem o modelo original não? Apesar de o intercepto não ter ficado muito correto, bem, existe um padrão nos dados, nos simulamos ele assim, e ele respondeu que existe esse padrão, e deu a forma mais ou menos correta.

Agora vamos ver para o segundo caso, onde a regressão linear não era de muita serventia.

> modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian) > summary(modelo_gam) Family: gaussian Link function: identity Formula: resposta_seno ~ s(preditor, k = 5, fx = T) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.9750 0.2667 18.66 3.46e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(preditor) 4 4 5.861 0.00181 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.401 Deviance explained = 48.4% GCV = 2.5598 Scale est. = 2.1332 n = 30

E de novo conseguimos detectar que algo está acontecendo, esses dados não apenas acaso puro, e vamos ver como é o modelo que ele está propondo gráficamente.

Agora sim não? Mas que magia é essa? Bem o GAM funciona a partir de regressões locais como se medíssemos como os dados mudam ao longo da nossa amostra. Algo como o modelos de Broken-Stick em regressão segmentada ou em piecewise.

Então a é algo como a gente pega uma janelinha nos dados e faz uma regressão.

Depois outra janelinha e faz outra regressão.

E a gente vai ficando com um monte dessas regressões.

Depois é só fazer uma função, tipo piecewise bem grandona e juntar todas elas.

No nosso modelo clássico de regressão, a gente tem

Y_i = \alpha + \beta \cdot X_i + \varepsilon_i

onde o \varepsilon_i=N(0,\sigma^2),

já que falamos que estamos um erro com distrubição normal dos resíduos.

A mudança que o GAM faz é

Y_i = \alpha + f (X_i) + \varepsilon_i

onde o \varepsilon_i=N(0,\sigma^2)

e o  f (X_i) é uma função que descreve os dados, e não mais um coeficiente como era na regressão.

So que você pode pensar, ei, então esse negocio vai ajustar qualquer coisa, qualquer conjunto de dados ele vai achar uma uma função que explica!

Não é tão simples assim, podemos fazer um teste, gerar dados sem padrão nenhum e ver o que ele encontra.

E para um conjunto de dados sem padrão.

Family: gaussian Link function: identity Formula: resposta_aleatoria ~ s(preditor, k = 5, fx = T) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.2277 0.1509 -1.509 0.144 Approximate significance of smooth terms: edf Ref.df F p-value s(preditor) 4 4 1.592 0.208 R-sq.(adj) = 0.0754 Deviance explained = 20.3% GCV = 0.81982 Scale est. = 0.68318 n = 30

Ele nos diz, não tem padrão nenhum aqui uai. O modelo simplesmente cobre tudo, não existe previsibilidade nenhuma.

Agora nem tudo são flores, e temos que tomar algumas decisões, sempre existe um preço a se pagar por algo tão legal assim, que é como definir a função que descreve os dados, veja que quando estavamo vendo como é feito ela, que é definida naquela função s(), separamos os dados em quadradinhos, fazendo pequenas regressões dentro deles certo, mas qual o tamanho ideal de quadradinho a usar? Isso é o que esta sendo definido naquele argumento k dentro da função s, e esta ligado ao nossos graus de liberdade, veja que quanto menor o quadradinho, menos dados por quadradinho então é mais facil definir errado a inclinação da reta, se tiver somente um ponto num quadradinho, nem tem como definir a inclinação da reta, agora se o quadradinho for muito grande, bem a regressão linear é como se fosse so um quadrado, então nao tem porque desse método, além de que podemos peder mudanças nos dados, então isso pode ser um problema, veja que no caso dos dados oscilando na função seno, se a gente usar um k muito baixo

modelo_gam<- gam(resposta_seno~s(preditor,k=2,fx=T),family = gaussian)

Avaliando os resíduos ainda vemos que algo não está correto.

Agora se aumentarmos o k
modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian)

E podemos ir aumentando até ficar aceitavel a distribuição dos residuos, a partir do momento que estiver aceitavel, aumentar mais é disperdicio pois estamos diminuindo os graus de confiança sem obter nenhum ganho em ter um modelo com melhor ajuste.

E é claro, alguém automatizou isso, então se você coloca -1 no k, a gente usa um algoritimo que define esse valor automaticamente.

Bem é isso ai, uma pequena introdução ao GAM, veja que aquela similaridade com a ideia de regressão não é so para ser bonitinho, isso permite trabalharmos com mais de um termo, interação entre termos, variáveis aleatória, bem tudo que da pra fazer com regressão da para fazer com GAM, sendo que regressão linear pode ser considerado um caso especial de GAM, que é quando a gente cai na função da reta. Eu acho que ando muito desanimado em casa, mas acho que preciso voltar a fazer posts com análise de dados que é muito legal pra ver se animo. O script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Zuur, AF 2012 Beginner’s Guide to Generalized Additive Models with R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
###Simulando Dados
set.seed(123)
preditor <- runif(30,0,4*pi)
resposta_pot<-5+2*(preditor)^2+rnorm(30,0,10)
resposta_seno<-5+2*sin(preditor)+rnorm(30)
 
##
jpeg("01.jpg")
plot(resposta_pot~preditor,frame=F,pch=19,xlab="Preditor", ylab="Resposta")
dev.off()
 
##
jpeg("02.jpg")
plot(resposta_pot~preditor,frame=F,pch=19,xlab="Preditor", ylab="Resposta")
modelo_linear <- lm(resposta_pot~preditor)
abline(modelo_linear, col="red")
dev.off()
 
##
summary(modelo_linear)
 
##
jpeg("03.jpg")
plot(resposta_pot~preditor,frame=F,pch=19,xlab="Preditor", ylab="Resposta")
curve(5+2*x^2,add=T,lwd=2,col="blue")
legend("topleft",col=c("blue","red"),lwd=c(2,1),legend=c("Realidade","Modelo"),bty="n")
 
modelo_linear <- lm(resposta_pot~preditor)
abline(modelo_linear, col="red")
novos_dados <- seq(0,4*pi,0.25)
pred_interval <- predict(modelo_linear, newdata=data.frame(preditor=novos_dados), interval="prediction",level = 0.95)
lines(novos_dados, pred_interval[,2], col="orange", lty=2)
lines(novos_dados, pred_interval[,3], col="orange", lty=2)
dev.off()
 
##
jpeg("04.jpg")
par(mfrow=c(2,2))
plot(modelo_linear)
dev.off()
 
##
jpeg("05.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12),xlab="Preditor",ylab="Resposta")
curve(5+2*sin(x),add=T,lwd=2,col="blue")
legend("topleft",col=c("blue","red"),lwd=c(2,1),legend=c("Realidade","Modelo"),bty="n")
 
modelo_linear <- lm(resposta_seno~preditor)
abline(modelo_linear, col="red")
novos_dados <- seq(0,4*pi,0.25)
pred_interval <- predict(modelo_linear, newdata=data.frame(preditor=novos_dados), interval="prediction",level = 0.95)
lines(novos_dados, pred_interval[,2], col="orange", lty=2)
lines(novos_dados, pred_interval[,3], col="orange", lty=2)
dev.off()
 
##
summary(modelo_linear)
 
##
jpeg("06.jpg")
par(mfrow=c(2,2))
plot(modelo_linear)
dev.off()
 
##Pacote utilizado  para o GAM
##install.packages("mgcv")
library(mgcv)
 
## Ajuste de modelo
modelo_gam<- gam(resposta_pot~s(preditor,k=5,fx=T),family = gaussian)
summary(modelo_gam)
 
## Visualização do modelo
jpeg("07.jpg")
plot(modelo_gam)
dev.off()
 
##
modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian)
summary(modelo_gam)
 
##
jpeg("08.jpg")
plot(modelo_gam)
dev.off()
 
##
jpeg("09.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12))
 
points(x=c(4,4,6,6,4),y=c(2,4,4,2,2),type="l",lty=2)
selecao <- preditor>=4 & preditor<=6
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],4,6,add=T,lwd=2,col="red")
dev.off()
 
##
jpeg("10.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12))
 
points(x=c(4,4,6,6,4)+3,y=c(2,4,4,2,2)*2+1,type="l",lty=2)
selecao <- preditor>=7 & preditor<=9
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],7,9,add=T,lwd=2,col="red")
dev.off()
 
##
jpeg("11.jpg")
plot(resposta_seno~preditor,frame=F,pch=19,ylim=c(0,12))
 
points(x=c(4,4,6,6,4),y=c(2,4,4,2,2),type="l",lty=2)
selecao <- preditor>=4 & preditor<=6
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],4,6,add=T,lwd=2,col="red")
 
points(x=c(4,4,6,6,4)+3,y=c(2,4,4,2,2)*2+1,type="l",lty=2)
selecao <- preditor>=7 & preditor<=9
modelo_linear <- lm(resposta_seno[selecao]~preditor[selecao])
curve(coef(modelo_linear)[1]+x*coef(modelo_linear)[2],7,9,add=T,lwd=2,col="red")
dev.off()
 
##Testando gam contra dados aleatórios
resposta_aleatoria <- rnorm(30)
 
##
jpeg("12.jpg")
plot(resposta_aleatoria~preditor,pch=19,frame=F)
dev.off()
 
modelo_gam<- gam(resposta_aleatoria~s(preditor,k=5,fx=T),family = gaussian)
summary(modelo_gam)
 
##
jpeg("13.jpg")
plot(modelo_gam)
dev.off()
 
##Definindo um valor de K
modelo_gam<- gam(resposta_seno~s(preditor,k=2,fx=T),family = gaussian)
jpeg("14.jpg")
hist(resid(modelo_gam))
dev.off()
 
##
modelo_gam<- gam(resposta_seno~s(preditor,k=5,fx=T),family = gaussian)
jpeg("15.jpg")
hist(resid(modelo_gam))
dev.off()
 
##Auto ajuste de k
modelo_gam<- gam(resposta_seno~s(preditor,k=-1,fx=T),family = gaussian)
jpeg("16.jpg")
hist(resid(modelo_gam))
dev.off()

Criando um gráfico com bandeiras dos países.

Esses tempos, um colega me pediu para fazer alguns gráficos tipo os que tem no google analytics, mais especificamente nesse estilo aqui.

É de uma versão mais antiga dessa figura, mas a ideia é essa, basicamente um gráfico de barras na vertical. A única parte difícil, foi as bandeiras, primeiro foi achar essas bandeiras, depois inserir na figura.

Para achar as figuras, a primeira coisa foi procurar algum pacote que tivesse isso pronto, achei que encontraria, mas não achei nenhum, então ao procurar pelas bandeiras, uma fonte seria baixar diretamente do wikipedia, que tem bandeiras para todos os paises, mas eu achei esse conjunto aqui, que inclusive está em um repositorio do github que podemos clonar, o que facilita bastante o trabalho, ali temos as bandeiras em formato png em vários tamanhos, eu preferi usar o de maior resolução.

então clonado o repositório com as bandeiras, eu fiz a figura com o seguinte código:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
library(png)
paises<-list.files("64")
 
##Figura country.png
 
country<-c("Brazil","United States","Portugal","Mozambique","Mexico","Germany","Peru","Australia","Cape Verde","Angola")
country_br<-c("Brasil","Estados Unidos","Portugal","Mozambique","México","Alemanha","Perú","Austrália","Cabo Verde","Angola")
sessions<-c(2113,259,85,63,54,49,38,37,36,36)
sessions_percentage<-c(60.39,7.4,2.43,1.8,1.54,1.4,1.09,1.06,1.03,1.03)
 
ampliacao <- 2
jpeg("country.jpg",width = ampliacao*480,height = ampliacao*480,pointsize = ampliacao*12)
escala<-barplot(rev(sessions_percentage),horiz=T,xlim=c(-150,100),xaxt="n",col="#76A7FA")
text(rep(0,10),escala,rev(sessions),pos=2)
text(rev(sessions_percentage)+14,escala,paste(rev(format(round(sessions_percentage, 2), nsmall = 2)),"%",sep=""))
text(rep(-130,10),escala,10:1,pos=2)
text(rep(-120,10),escala,rev(country_br),pos=4)
mtext(c("País","Sessões","% Sessões"),3,0,at=c(-125,-15,35))
 
list.files()
 
for(i in 1:10){
    country[i]
    posicao<-agrep(country[i],paises)
    paises[posicao]
    bandeira<-readPNG(paste("64/",paises[posicao],sep=""))
    rasterImage(bandeira,-133,rev(escala)[i]-0.5,-117,rev(escala)[i]+0.5)
}
dev.off()

Que gera a seguinte figura.

Bem legal eu achei, agora vou fazer alguns comentários do meu código.

Primeiro:

1
2
country<-c("Brazil","United States","Portugal","Mozambique","Mexico","Germany","Peru","Australia","Cape Verde","Angola")
country_br<-c("Brasil","Estados Unidos","Portugal","Mozambique","México","Alemanha","Perú","Austrália","Cabo Verde","Angola")

Bem como são poucos dados, eu digitei o nomes do países em inglês e em português, isso porque o nome em inglês eu uso para achar o nome do arquivo da bandeira, que é o nome do pais, e em português para usar na figura em si.

Depois disso para facilitar a mudança do tamanho, eu costumo guardar uma variável para o tamanho, que multiplica a largura, a altura e o tamanho do ponto, assim só altero ela para mudar o tamanho da figura, a sua resolução.

1
2
3
4
5
6
ampliacao <- 2
jpeg("country.jpg",width = ampliacao*480,height = ampliacao*480,pointsize = ampliacao*12)
 
##Figura
 
dev.off()

A figura em si é um barplot na vertical, a única coisa é que mudo bastante as margens, o tamanho do eixos de forma a caber as outras informações e bandeiras, que são colocadas com a função text para dentro da figura e mtext para a parte de fora. Depois disso, talvez essa parte ainda tenha algo de interesse.

1
2
3
4
5
6
7
for(i in 1:10){
    country[i]
    posicao<-agrep(country[i],paises)
    paises[posicao]
    bandeira<-readPNG(paste("64/",paises[posicao],sep=""))
    rasterImage(bandeira,-133,rev(escala)[i]-0.5,-117,rev(escala)[i]+0.5)
}

Então eu pego o nome do pais em inglês, e uso agrep na lista de nomes de arquivos, agrep, porque é um grep, ou seja acha a palavra na lista, mas permite alguns erros de caracteres, que da para regular com argumentos, mas então com o agrep a gente vai achar o nome do arquivo, com ele a gente le a figura png com o readPNG, que le o raster da figura, então com rasterImage a gente adiciona ela a figura, aqui a gente tem que colocar os quatro pontos da imagem para adicionar ela, veja se da até para colocar a imagem espelhada, feito isso ta pronto. Para economizar, eu fiz um loop para todas as bandeiras.

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail. Um post bem simples, mas é o tipo de código que é legal guardar, pois pode ser útil copiar e colar um dia.

Variância e a desigualdade de Chebyshev

Como falamos antes de média aqui, vamos falar de Variância agora.

Variância é definida como:

Var(X)=E[(X-\mu)^2]

Sendo que esse E é o valor esperado, que é

 E[X]= \sum\limits_{i=1}^\infty x_i \cdot p(x_i)

e o \mu é o centro da distribuição, a média, como vimos aqui

E essa é a forma geral de calcular a variância, que é o espalhamento dos dados em torno da média.

Uma forma mais fácil de calcular é

Var(X)=E[X^2]-E[X]^2

A raiz quadrada da variância é o famoso desvio padrão, que como não está em unidades quadradas, vai estar na mesma unidade de medida dos dados originais. Veja que se a gente tem uma constante multiplicando X

Var(aX)= a^2 \cdot Var(X)

Ela sai da variância ao quadrado.

Para entender melhor vamos calcular a variância de uma jogada de dados, o dado comum de 6 faces.

Sabemos que a valor esperado é:

E[X]=3.5

e

E[X^2]=1^2\cdot\frac{1}{6}+2^2\cdot\frac{1}{6}+3^2\cdot\frac{1}{6}+4^2\cdot\frac{1}{6}+5^2\cdot\frac{1}{6}+6^2\cdot\frac{1}{6}=15.17

Então para a variância temos:

Var(X)=E[X^2]-E[X]^2 = 15.17-3.5^2=2.92

Veja que se a gente usar:

Var(X)=E[(X-\mu)^2]= (1-3.5)^2\cdot\frac{1}{6}+(2-3.5)^2\cdot\frac{1}{6}+(3-3.5)^2\cdot\frac{1}{6}+(4-3.5)^2\cdot\frac{1}{6}+(5-3.5)^2\cdot\frac{1}{6}+(6-3.5)^2\cdot\frac{1}{6}

Como a média está em todas as multiplicações, podemos fatorar ela, que é da onde vem a fórmula anterior, lembrando que a média está dentro de um quadrado, então não da para ir tirando ela diretamente.

Agora qual a variância de uma jogada de moeda com uma chance p para cara.

E[X]=0 \cdot (1-p) + 1 \cdot p = p


E[X^2]= E[X] = p


Var(X)=E[X^2]-E[X]^2=p-p^2=p(1-p)

E essa fórmula é bem conhecida, mesmo sem saber da onde ela vinha, e veja que ela é maximizada no p=0.5

Ou seja, a maior entropia, menos previsibilidade, está numa moeda honesta.

Então suponha que uma variável aleatória X em que O \leq X \leq 1 e E[X]=p

Veja que X^2 \leq X pode ser igual quando pensamos em 1 por exemplo, já que um ao quadrado da 1, mas fora isso, sempre X^2 vai ser maior, então E[X^2]\leq E[X]=p

Assim, Var(X) = E[X^2]-E[X]^2 \leq E[X]-E[X]^2 = p(1-p)

Dessa forma, a variância de Bernoulli é a maior possível para uma variável aleatória entre 0 e 1. O que nos leva a inequação de Chebyshev, que é útil para interpretar de forma geral variâncias.

A inequação mostra que

P( X-\mu \leq k\sigma) \leq \frac{1}{k^2}

Por exemplo, a probabilidade que uma variável aleatória esteja além de k desvios padrões é menos de \frac{1}{k^2}

 2 \sigma = 25\%


 3 \sigma = 11\%


 4 \sigma = 6\%

Mas veja que esse é apenas um limite geral, a probabilidade pode ser bem menor, mas esse é um limite superior independente da distribuição, e essa generalização que torna essa inequação importante, principalmente para afirmações de forma geral.

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e esses post está a uns dias parados, então resolvi finalizar logo. Além disso tem uma prova para essa inequação, mas achei a prova muito difícil.

Referência:
Coursera, curso Mathematical Biostatistics Boot Camp 1 do Brian Caffo

Valor esperado (expected value), a famosa média.

O valor esperado (expected value), esperança matemática ou média de uma variável aleatória é o centro da sua distribuição.

Para uma variável aleatória discreta X com uma PMF p(x), ela pode ser definida da seguinte forma.

E[X] = \sum \limits_x x \cdot p(x)

Ou seja, cada valor multiplicado por sua probabilidade.
E[X] representa o centro de massa de uma coleção de locais e seus pesos \{ x,p(x)  \}

Isso é meio idiota, mas as vezes ao esquecer essas coisas simples, agente não entende coisas mais complexas, mais para frente. Além disso da para ter uma intuição aqui, quando falamos que média ou valor esperado é o centro de massa, é porque ele é o ponto onde podemos equilibrar a distribuição.

A partir da formula ali em cima, podemos fazer uma função no R para calcular então o valor esperado:

1
2
3
4
5
6
7
esperanca<-function(vetor){
    esperanca<-0
    for(i in 1:length(vetor)){
        esperanca<-esperanca+ i * (vetor[i]/sum(vetor))
    }
    return(esperanca)
}

Veja que é o somatório, fazendo a multiplicação, não é a forma mais rápida de implementar, mas eu acho a leitura dos loops dessa forma mais simples do que usando operações vetoriais, mas tudo bem, suponha que temos uma distribuição com 7 possibilidades de valores X=\{1,2,3,4,5,6,7 \}, vamos ver como fica o valor esperado para diferentes probabilidades para cada caso.

Estamos tentando equilibrar a distribuição na ponta daquela seta, veja nos exemplos o que acontece conforme mudamos as probabilidades, principalmente, veja o quanto no caso C, uma probabilidade alta puxa para um lado o valor esperado, mas veja que uma probabilidade baixa, como no caso D, não puxa muito a seta para o lado dela.

Agora muitas vezes o que a gente vê as pessoas usando, é essa fórmula, mas com o quadrado do X.

E[X^2] = \sum \limits_x x^2 \cdot p(x)

Isso resolve o problemas dos valores negativos, ja que todo valor ao quadrado será positivo, mas mais pra frente veremos melhor isso.

Então a aplicação da fórmula é direta, para um exemplo simples, podemos pensar no lançar de uma moeda honesta, onde X= \{ 0,1 \} e p(0)=0.5 e p(1)=0.5, então

E[X] = \sum \limits_x x \cdot p(x) =  0 \cdot 0.5 + 1 \cdot 0.5 = 0.5

O valor esperado para uma moeda honesta é 0.5, veja que é interessante pensar, que 0.5 não é um valor válido para a jogada da moeda, mas é o valor esperado, ou seja o valor esperado não é necessariamente o valor que mais acontece, e pode nem ser um valor que acontece, como nesse caso.

Para uma variável contínua, com densidade f, o valor esperado vai ser definido da seguinte forma.

E[X]=\int_{-\infty}^{\infty} t \cdot f(t)dt

Que vem da definição da física de centro de massa, e lembrando que mesmo que a distribuição só ocupe uma parte dos valores reais, o resto da distribuição terá probabilidade zero.

Vamos ver o caso da distribuição uniforme de mínimo 0 e máximo 1.

Primeiro, não se engane com o gráfico, veja que temos um quadrado, com o lado 1 e topo 1, então calculando a integral, ou a área desse quadrado, é bem intuitivo que ele é 1, então ele tem área um, e nenhum valor menor que zero, então é uma distribuição válida, é um PMF válido e o valor esperado vai ser:

E[X]=\int_{0}^{1} x \cdot dx = \frac{x^2}{2} \Big|_0^1 = \frac{1}{2}

Que é o quantile de 50%, ou 0.5

> qunif(0.5,min=0,max=1) [1] 0.5

Agora, o legal é que existem algumas regras, que podem ajudar bastantes, regras quanto ao valor esperado, que é um operador linear.
Se a e b não são valores fixos e X e Y são duas variáveis aleatórias, temos que:

E[a \cdot X+b]=a \cdot E[X]+b

e

 E[X+Y] = E[X] + E[Y]

Em geral, se g é uma função não linear,

 E[g(x)] \neq g(E[X])

Por exemplo E[X^2] \neq E[X]^2

Vamos supor que vamos lançar dois dados honestos, qual o valor esperado da soma dos resultado?

[E[Dado_1+Dado_2]=E[Dado_1]+E[Dado_2]]
Sabemos que um dado tem valor esperado de 3.5 pois:

E[Dado]=\frac{1}{6} \cdot 1 +\frac{1}{6} \cdot 2 +\frac{1}{6} \cdot 3 +\frac{1}{6} \cdot 4 +\frac{1}{6} \cdot 5 +\frac{1}{6} \cdot 6 = 3.5

logo temos que:

[E[Dado_1+Dado_2]=E[Dado_1]+E[Dado_2]]=3.5+3.5=7

Certo, e o valor esperado da média da jogada de dois dados?

[E[\frac{Dado_1+Dado_2}{2}]=\frac{1}{2} \cdot(E[Dado_1]+E[Dado_2]])=\frac{1}{2} \cdot(3.5+3.5)=3.5

Que interessante, o valor esperado da média de dois dados é igual ao valor esperado de um dado, e isso é válido para qualquer número de dados, ou qualquer coleção de variáveis.

Seja X_i parai=1,2,\dots,n uma coleçao de variáveis aleatórias, cada uma de uma distribuição com média \mu, o valor esperado da amostra média de X_i será:

E[\frac{1}{n} \sum\limits_{i=1}^n X_i ]

\frac{1}{n} \cdot E[ \sum\limits_{i=1}^n X_i ]

\frac{1}{n} \cdot  \sum\limits_{i=1}^n E[X_i ]

E como sabemos que E[X_i ]=\mu

\frac{1}{n} \cdot  \sum\limits_{i=1}^n \mu = \mu

Assim, o valor esperado da média da amostra é a média da população que estamos tentando estimar. Quando o valor esperado de um estimador é o que estamos tentando estimar, nós dizemos que esse estimador é não enviesado, que é a ideia que seguimos na estatística frequentista certo, podemos trabalhar com amostras, e estimar coisas dela, porque a estimativa do valor esperado da média deve ser o valor estimado da população inteira, e desse monte de resultado talvez pouco intuitivos, temos um resultado forte aqui.

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Coursera, curso Mathematical Biostatistics Boot Camp 1 do Brian Caffo

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
##Função para calcular o valor esperado
esperanca<-function(vetor){
    esperanca<-0
    for(i in 1:length(vetor)){
        esperanca<-esperanca+ i * (vetor[i]/sum(vetor))
    }
    return(esperanca)
}
 
##Figura 1
par(mfrow=c(2,2))
valores<-c(4,4,0,0,0,4,4)
plot(1:length(valores),valores,ylim=c(0,10),ylab="Probabilidade",yaxt="n",frame=F,xaxt="n",xlab="",pch=19,cex=2)
for(i in 1:length(valores)){
    lines(c(i,i),c(0,valores[i]),lwd=4)
}
axis(2,at=seq(0,10,2),labels=round(seq(0,10,2)/16,2),las=1)
axis(1,at=1:length(valores),labels=NA)
valor_esperado <- esperanca(valores)
arrows(valor_esperado, -2,valor_esperado,-1, xpd = TRUE,lwd=10)
legend("top",legend="A",bty="n",cex =2)
 
valores<-c(4,0,0,0,4,4,4)
plot(1:length(valores),valores,ylim=c(0,10),ylab="Probabilidade",yaxt="n",frame=F,xaxt="n",xlab="",pch=19,cex=2)
for(i in 1:length(valores)){
    lines(c(i,i),c(0,valores[i]),lwd=4)
}
axis(2,at=seq(0,10,2),labels=round(seq(0,10,2)/16,2),las=1)
axis(1,at=1:length(valores),labels=NA)
valor_esperado <- esperanca(valores)
arrows(valor_esperado, -2,valor_esperado,-1, xpd = TRUE,lwd=10)
legend("top",legend="B",bty="n",cex =2)
 
valores<-c(9,0,0,0,4,2,1)
plot(1:length(valores),valores,ylim=c(0,10),ylab="Probabilidade",yaxt="n",frame=F,xaxt="n",xlab="",pch=19,cex=2)
for(i in 1:length(valores)){
    lines(c(i,i),c(0,valores[i]),lwd=4)
}
axis(2,at=seq(0,10,2),labels=round(seq(0,10,2)/16,2),las=1)
axis(1,at=1:length(valores),labels=NA)
valor_esperado <- esperanca(valores)
arrows(valor_esperado, -2,valor_esperado,-1, xpd = TRUE,lwd=10)
legend("top",legend="C",bty="n",cex =2)
 
valores<-c(1,0,0,0,1,5,10)
plot(1:length(valores),valores,ylim=c(0,10),ylab="Probabilidade",yaxt="n",frame=F,xaxt="n",xlab="",pch=19,cex=2)
for(i in 1:length(valores)){
    lines(c(i,i),c(0,valores[i]),lwd=4)
}
axis(2,at=seq(0,10,2),labels=round(seq(0,10,2)/16,2),las=1)
axis(1,at=1:length(valores),labels=NA)
valor_esperado <- esperanca(valores)
arrows(valor_esperado, -2,valor_esperado,-1, xpd = TRUE,lwd=10)
legend("top",legend="D",bty="n",cex =2)
 
##Figura 2
valores<-c(0.5,0.5)
plot(1:length(valores),valores,ylim=c(0,1),ylab="Probabilidade",yaxt="n",frame=F,xaxt="n",xlab="",pch=19,cex=2,main="Moeda {0,1}")
for(i in 1:length(valores)){
    lines(c(i,i),c(0,valores[i]),lwd=4)
}
axis(2,at=seq(0,1,0.1),las=1)
axis(1,at=1:length(valores),labels=NA)
valor_esperado <- esperanca(valores)
arrows(valor_esperado, -2,valor_esperado,-0.05, xpd = TRUE,lwd=10)
 
##Figura 3
curve(dunif(x),-0.5,1.5,frame=F,xlab="",ylab="Densidade",main="Distribuição uniforme")
 
##Valor esperado da distribuição uniforme, quantile de 50%
qunif(0.5,min=0,max=1)

CDFs e Quantiles

Bem, então distribuições estatísticas são funções matemáticas que definem a probabilidade de certos eventos ocorrerem, um modelo para o comportamento de variáveis aleatórias. No post anterior vimos como são as PDF e PMF, ou funções de densidade, que nos mostram como essas probabilidades são, e para a distribuição exponencial vimos como com uma integral pode ser usada para definir a probabilidade de um determinado evento.

Essa funções para calcular essas areas de interesses são tão importantes que tem nomes especiais, que são as cumulative distribution function (CDF) e survival function.

Uma CDF de uma variável aleatória X é definida como a função

F(x)=P(X \leq x)

Essa definição se aplica independente de X ser discreto ou contínuo.
Já a survival function de uma variável X é definida como
S(x)=P(X>x)
Veja que como já definimos F(x), a survival function será S(x)=1-F(x), porque S(x)+F(x)=1, você se lembra como toda a área da probabilidade tem que ser 100% ou 1, então é isso.
Além disso, para variáveis aleatórias contínuas, o PDF é a derivada do CDF.

Vamos tentar ver isso baseado no exemplo do post passado aqui onde x=6.

Nosso PDF era
  f(x) = \left\{  \begin{array}{ll}  \frac{e^{-x/5}}{5} & \quad x\ \textgreater \ 0 \\  0 & \quad x \leq 0  \end{array}  \right.

Então a CDF, e a survival function são

Assim a surival function vai ser

\int_{x}^{\infty} \frac{e^{(-t/5)}}{5} dt = -e^{(-t/5)} \Big|_x^\infty=-e^{(-x/5)}

Assim sabemos que

F(x)=1-S(x)=1-e^{(-x/5)}

Porém podemos notar que podemos voltar na PDF com a derivada

f(x)=F'(x)=\frac{d}{dx}(1-e^{(-x/5)})=\frac{e^{(-x/5)}}{5}

Mas calcular isso no R, basta usar o pexp

> pexp(6,1/5,lower.tail = F) [1] 0.3011942 > 1-pexp(6,1/5,lower.tail = F) [1] 0.6988058

E podemos colocar esses valores na figura

Então a área vermelha é 30% da probabilidade e a área branca é 70% da probabilidade.

Assim, eu falo um ponto, em anos, e eu te digo quanto é a probabilidade da pessoa estar viva ou morrer, segundo nosso modelo. Agora suponha que a minha pergunta fosse, com quantos anos, após o exame, a pessoa acumula 50% de chance de morrer? Para essa pergunta a gente precisa do quantile.

O \alpha^{th} quantile de uma distribuição com função de distribuição F é o ponto x_\alpha tal que

F(x_\alpha)=\alpha

Um percentile é simplismente um quatile \alpha representado por porcentagem.

A mediana é o 50^{th} percentile

Para achar esse quantile, o que queremos é

0.5=F(x)

0.5=1-e^{(-x/5)}

x=-log(0.5) \cdot 5 \approx  3.47

Essa função esta implementada com o q antes do nome das distribuições, então nesse caso é qexp para achar o 50%

> qexp(0.5,1/5) [1] 3.465736

E para entender esse número, podemos olhar ele na nossa figura

Então quando a gente olha ?Distributions no R, vemos que cada distribuição tem 4 funções básicas,

  • A PDF, que começa com d
  • A CDF, que começa com p
  • O quantile que começa com q
  • E um gerador de números aleatórios para cada distribuição, que começa com r

E agora a gente sabe o que cada uma dessas funçõezinhas fazem 🙂

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Coursera, curso Mathematical Biostatistics Boot Camp 1 do Brian Caffo

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
 
##jpeg("01.jpg")
curve(dexp(x,1/5),0,20,frame=F,xlab="Tempo em anos",ylab="Probabilidade")
polygon(c(seq(6,20,0.1),seq(20,6,-0.1)),c(dexp(seq(6,20,0.1),1/5),rep(0,141)),col="red")
text(8,0.02,"S(x)")
text(4,0.02,"F(x)")
##dev.off()
 
pexp(6,1/5,lower.tail = F)
 
1-pexp(6,1/5,lower.tail = F)
 
##jpeg("02.jpg")
curve(dexp(x,1/5),0,20,frame=F,xlab="Tempo em anos",ylab="Probabilidade")
polygon(c(seq(6,20,0.1),seq(20,6,-0.1)),c(dexp(seq(6,20,0.1),1/5),rep(0,141)),col="red")
text(8,0.02,paste("S(x)=",round(pexp(6,1/5,lower.tail = F),3),sep=""))
text(4,0.02,paste("F(x)=",1-round(pexp(6,1/5,lower.tail = F),3),sep=""))
##dev.off()
 
qexp(0.5,1/5)
 
##jpeg("03.jpg")
q50<-qexp(0.5,1/5)
curve(dexp(x,1/5),0,20,frame=F,xlab="Tempo em anos",ylab="Probabilidade")
polygon(c(seq(q50,20,0.1),seq(20,q50,-0.1)),c(dexp(seq(q50,20,0.1),1/5),rep(0,166)),col="red")
text(q50-2,0.02,paste("S(x)=",round(pexp(q50,1/5,lower.tail = F),3),sep=""))
text(q50+2,0.02,paste("F(x)=",1-round(pexp(q50,1/5,lower.tail = F),3),sep=""))
mtext(round(q50,2),side=1,at=q50)
abline(v=q50,col="blue",lty=3,lwd=2)
##dev.off()