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()

Teorema de Bayes e por que não devemos discutir com com quem tem certeza absoluta de alguma coisa.

Eu acho o teorema de Bayes muito legal, ja falamos dele muitas vezes, como aqui e aqui, além do que está espalhado em muitos outros posts, mas é sempre possível aprender algo novo.

O teorema de bayes tenta de dar a probabilidade de você ter uma doença (H, de hypothesis) dado alguma evidência, que veio de um evento (E, de event), como receber um teste positivo de um exame para detectar alguma doença, a probabilidade de H dado que E aconteceu p(H|E).

Para isso você precisa saber, qual o prior, a chance de você ter a doença p(H) e multiplicar pela probabilidade do evento acontecer dado que a hipótese é verdadeira p(E|H), e então dividir pela probabilidade total do evento ocorrer, p(E), a probabilidade de testar positivo.

p(H|E)=\frac{ p(E|H) \cdot p(H)}{p(E)}

Mas o que é p(E), bem é uma combinação de ter a doença e e ser corretamente identificado por ela p(H) \cdot P(E|H) mais a chance de não ter a doença e ser falsamente identificado. p(-H) \cdot P(E|-H), ou seja

p(H|E)=\frac{ p(E|H) \cdot p(H)}{p(H) \cdot P(E|H)+p(-H) \cdot P(E|-H)}

Agora, o exemplo clássico do exame, sabemos que você faz um exame, o resultado volta e você tem um positivo para uma rara doença que acomete a população, somente 0.1% da população tem essa doença, então você se pergunta o quão certo eu posso ser disso, bem o teste identifica corretamente a doença 99% das vezes, só errando 1%, então é 99% certo que você tem a doença? Vamos usar o teorema de Bayes.

p(H)=0.001

p(E|H)=0.99

P(E|-H)=1-p(E|H)=0.01

O que da.
p(H|E)=\frac{ 0.99 \cdot 0.001}{0.001 \cdot 0.99 + (1-0.001) \cdot 0.01}=9%

E parece perdido esse resultado, o exame é 99% certeiro para identificar a doença, ainda assim você tem um positivo e agora acha que tem 9% de chance de ter a doença?Visualizando fica mais fácil.

Primeiro, sabemos que a doença é rara, ou seja, somente 0.1% da população tem ela, isso é uma pessoa em mil.

Ou seja, se os 1000 pontinhos são a população, somente uma, a bolinha vermelha é doente nessa população. Mas o teste é eficaz em 99% das vezes, ele erra 1% das vezes, errar é dizer que uma pessoa sadia esta doente aqui. Então se testarmos as outras 999 pessoas que não são doentes, o teste vai falhar em 1% dos casos, ou seja, 0.01*999, que da 9.99, ou seja, 10 pessoas vão receber positivo, mesmo sendo sadias.

Ou seja, p(H|E), agora que tiramos as bolinhas pretas da jogada, qual a chance de eu estar doente?

O que da aproximadamente 9%, uma bolinha doente num mar de 11 bolinhas, 1/11.

> 1/11 [1] 0.09090909

Basicamente o mesmo resultado de aplicar o teorema de Bayes

> (0.99*prior)/(prior*0.99+(1-prior)*0.01) [1] 0.09016393

Agora o por que não devemos discutir com com quem tem certeza absoluta de alguma coisa?

Suponha que você acha que essa doença é crendice popular, que não tem nada a ver, é da cabeça das pessoas, você acha que seu prior então é zero.

> prior<-0.0 > (0.99*prior)/(prior*0.99+(1-prior)*0.01) [1] 0

Veja que no final você continua com a mesma probabilidade posterior, se a doença não existe, não tem como alguém ter ela.

Agora suponha que você acha que existe uma conspiração do governo, e que na verdade todo mundo ta doente, e o governo nos engana, então o prior é 100%, se usarmos o teorema de Bayes temos:

> prior<-1.0 > (0.99*prior)/(prior*0.99+(1-prior)*0.01) [1] 1

Ou seja, todo mundo já ta doente, não precisa de evidência mostrando que as pessoas estão doentes, nos dois extremos, 0% ou 100% de prior, nossa opinião não muda nada, por isso que não devemos atribuir 0% nem 100% de chance para nada, porque se ja temos essa certeza, evidência de algo contrário ou não ao que achamos é inútil, não vai mudar nossa perspectiva.

Bem é isso ai, tentando voltar a atividade, porque fiquei umas semanas desanimado com a vida, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Video: The Bayesian Trap do Veritasium.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
set.seed(123)
coordenadas<-expand.grid(1:50,1:20)
pessoas<-rep("black",1000)
pessoas[sample(1:1000,1)]<-"red"
 
plot(coordenadas[,1],coordenadas[,2],frame=F,xaxt="n",yaxt="n",xlab="1000 Pessoas",ylab="",pch=19,cex=1.2,col=pessoas)
 
 
pessoas[sample(1:1000,10)]<-"blue"
plot(coordenadas[,1],coordenadas[,2],frame=F,xaxt="n",yaxt="n",xlab="1000 Pessoas",ylab="",pch=19,cex=1.2,col=pessoas)
 
pessoas[pessoas=="black"]<-"white"
plot(coordenadas[,1],coordenadas[,2],frame=F,xaxt="n",yaxt="n",xlab="11 Pessoas",ylab="",pch=19,cex=1.2,col=pessoas)
 
prior<-0.001
(0.99*prior)/(prior*0.99+(1-prior)*0.01)
 
1/11
 
prior<-0.0
(0.99*prior)/(prior*0.99+(1-prior)*0.01)
 
prior<-1.0
(0.99*prior)/(prior*0.99+(1-prior)*0.01)

PMFs e PDFs

Sempre estamos falando em distribuições estatísticas, mas nunca vamos muito a fundo nelas, sempre paramos no meio do caminho, mas talvez tentar entender um pouco mais sobre distribuições seja algo valido. Podemos começar com Probability mass function (PMF) e Probability density function (PDF), que são funções matemáticas que servem para mapear regras de probabilidades de variáveis aleatórias. No caso PMF para variáveis discretas e PDF para variáveis contínuas.

Então uma função de massa de probabilidade (em português), ou PMF (como a maioria da literatura é em inglês, eu acho melhor usar os termos em inglês mesmo, gera menos confusão com livros) é uma função que você pluga um número nela e temos a probabilidade correspondente da variável aleatória do valor avaliado na função. Mas para ser uma função PMF valida, existem algumas propriedades que a função deve satisfazer. Considerando x como o conjunto de todas as possibilidades, como por exemplo para uma moeda x=\{cara,coroa \} ou um dado em que x=\{1,2,3,4,5,6 \}, temos que:

\forall x : p(x) \ge 0

Ou seja, para todo x, a probabilidade tem que ser zero ou maior que zero, porque não existe probabilidade negativa.

\sum_x p(x) = 1

E bem, se x é o conjunto de todas as possibilidades, a probabilidade de todas as possibilidades que que ser 100%, no caso 1.

Então, para o caso da moeda x=\{cara,coroa \}, pensando em cara como 0 e coroa como 1 x=\{0,1 \}, podemos construir a PMF da assim:

p(x)= \frac{1}{2}^x \cdot \frac{1}{2}^{(1-x)}

Veja que se sair cara, cara é igual a 0, da

p(0)= \frac{1}{2}^0 \cdot \frac{1}{2}^{(1-0)} p(0)= 1 \cdot \frac{1}{2}^{1} p(0)= \frac{1}{2}

Agora para coroa, que é 1, teremos o seguinte
p(1)= \frac{1}{2}^1 \cdot \frac{1}{2}^{(1-1)}

p(1)= \frac{1}{2} \cdot \frac{1}{2}^{(0)} p(1)= \frac{1}{2} \cdot 1 p(1)= \frac{1}{2}

E se somamos as duas probabilidades, de cara e coroa da 1, parece bem idiota, mas é o tipo de coisa idiota que nunca é ensinado eu acho.
Esse é o caso para uma moeda honesta, mas de qualquer forma, uma moeda, ou qualquer variável com probabilidade diferente, podemos generalizar da seguinte forma.

p(x)= (\theta)^x \cdot (1-\theta)^{(1-x)}

Onde \theta é a probabilidade dar alguns dos resultados, e temos um PMF para uma moeda não honesta, que não tem 50% de chance para cara ou coroa, e podemos ver que ela tem a soma de 1, porque:

Para temos x=\{0,1 \}
então
\sum_x p(x) = 1
vai ser
p(0) + p(1) = 1

(\theta)^0 \cdot (1-\theta)^{(1-0)} + (\theta)^1 \cdot (1-\theta)^{(1-1)} = 1

 1-\theta + \theta = 1

 1 = 1

Certo, mas isso foi para variáveis discretas, no caso de variáveis contínuas, como medidas de peso, altura, quantidade de energia de um individuo gasta temos as PDF, nesse caso, a área sobre a PDF corresponde a probabilidade de interesse, sendo que talvez a PDF mais famosa seja a da distribuição normal.

Para uma função ser uma PDF, ela tem que satisfazer a condições similares a da PMF, nesse caso

\forall x : f(x) \ge 0

Nada de probabilidade negativa, isso não faria sentido e

\int_{-\infty}^{\infty} f(x)dx = 1

Veja que aqui, temos que a área é 1, mas definimos a função sobre toda a reta dos números reais, porque mesmo que a função só cubra uma parte dos reais, sei la de 0 a 10, basta estar no zero, para outras posições.

Vamos ver um exemplo, suponha que queremos representar o tempo em anos do diagnostico até a morte de pessoas com um tipo de câncer específico que assume uma função de densidade da seguinte forma:

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

Para facilitar a vida, vamos escrever ela no R e desenhar que fica mais fácil visualizar.

1
2
3
fx <- function(x) {
    return(ifelse(x>0,exp(-x/5)/5,0))
}

E usando curve para plotar temos.

Ou seja, a chance de você morrer a alguns anos dos diagnostico é alta, e viver muito tempo, e morrer depois de 20 anos do diagnostico é baixo, bem intuitivo certo, porque cancer é uma doença complicada. Agora veja que o tempo negativo não tem probabilidade, porque se a pessoa não está viva para o diagnostico, não existe o tempo em anos do diagnostico até a morte.

Bem da para ver que essa é uma pdf válida porque e elevado a qualquer valor é sempre positivo, e

\int_{0}^{\infty} f(x)dx

\int_{0}^{\infty} \frac{e^{(-x/5)}}{5} dx

 \frac{1}{5}\int_{0}^{\infty} e^{(-x/5)} dx

-e^{(-x/5)} \Big|_0^\infty

Mas eu não sei fazer integral direito, mas podemos integrar nossa função no R, para ter certeza

> integrate(fx,0,Inf) 1 with absolute error < 2e-07

Ta mas como usa essa pdf, como aqui os valores são continuos, temos que fazer a pergunta certa, por exemplo, qual a chance de uma pessoa sobreviver mais de 6 anos.

Então a gente quer saber p(x \geq 6), que da.

p(x \geq 6)=\int_{0}^{\infty} \frac{e^{(-t/5)}}{5} dt

p(x \geq 6)=-e^{(-t/5)} \Big|_6^\infty

p(x \geq 6)=-e^{(-6/5)} \approx 0.301

A gente pode conferir com o integrate do R

> integrate(fx,6,Inf) 0.3011942 with absolute error < 3.8e-05

Ou, melhor ainda, a gente pode começar a ver que isso ta tudo pronto no R, sendo que essa é a distribuição exponencial, então a gente calcula com pexp

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

Sendo que no caso do pexp colocamos lower.tail = F, porque no nosso caso, quando x é menor que zero, temos zero, so queremos a parte positiva da cauda da distribuição.

Se você digitar no R

1
?Distributions

Você vai ver todas as distribuições implementadas no R por padrão, e todas tem a suas respectivas PDFs, dexp para a exponencial, dnorm para a normal, dgamma para a gamma, e agora a gente entende melhor o que são elas.

Então PMF e PDF são modelos que descrevem uma variável aleatória, que como não podemos saber seu valor, podemos apenas atribuir uma probabilidade

Bem é isso ai, hoje não precisa de script certo, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e vejam o curso da onde eu tirei essas coisas na referência, muito legal.

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
##Definindo uma pdf
fx <- function(x) {
    return(ifelse(x>0,exp(-x/5)/5,0))
}
 
##Gráfico da função
curve(fx(x),-20,20,frame=F,xlab="Tempo em anos",ylab="Probabilidade")
 
##integral de fx
integrate(fx,0,Inf)
 
integrate(fx,6,Inf)
 
pexp(6,1/5,lower.tail = F)

Regressão segmentada ou piecewise.

Depois de um tempo sumido, espero voltar finalmente a postar aqui. Estou devagar para estudar ultimamente, é de dezembro meu último post aqui, mas agora que acabou minha faculdade, vou voltar a estudar novas coisas e postar regularmente, existe muito coisa que pensei em postar e parei.

Mas para retomar as atividades, vamos falar sobre a piecewise regression, ou regressão segmentada. A gente ta acostumado a modelos de diferenças de médias, tipo anova, e modelos de regressão linear, ou a combinação desses dois. Eles resolvem bem quase todos os casos, mas as vezes a gente já sabe um pouco mais sobre o sistema que estamos estudando, e podemos fazer melhor, utilizando esse conhecimento.

Por exemplo, sabemos que plantas muitas vezes são limitadas, no seu crescimento, pelo mínimo de alguns nutrientes, a lei de Liebig. Podemos expandir esse conceito e pensar assim, enquanto não há um mínimo de nutrientes no solo, o crescimento é estagnado, mas a partir de algum ponto começamos a ter um incremento no crescimento de uma espécie dado um acréscimo nos nutrientes.

Podemos simular alguns dados para essa situação, então imaginando que temos uma medida dos nutrientes num vaso e do tamanho da planta após um determinado tempo fixo, podemos visualizar esses dados da seguinte forma.

Veja então que no eixo x temos a concentração de nutrientes, e no eixo y temos o tamanho das plantulas, então enquanto não temos uma certa concentração parece que todo mundo tem o mesmo tamanho, mas a partir de uma concentração, quanto maior a concentração, maior ficam as plantulas. Se a gente pensar, isso pode ocorrer em um bocado de situações, o que interessa é que, a partir de um ponto, as coisas mudam, algo começa a ocorrer, ou para, aqui temos um ponto crítico, algo entre o 4 e o 6 em que abaixo dele, a concentração de nutrientes é irrelevante.

Podemos representar isso matematicamente como uma equação do tipo piecewise, nesse caso, o que usamos para simular os dados foi:

  f(x) = \left\{          \begin{array}{ll}              8 & \quad x\ \textless \  5 \\              1+2x & \quad x \geq 5          \end{array}      \right.

Sendo x o nosso preditor, concentração de nutriente e f(x) o tamanho, para a parte determinística do modelo.
Esses são os valores que usamos para gerar essas dados de exemplo, que no caso é partido no 5, antes do 5 temos uma reta f(x)=8+0x(eu omiti o 0x, porque tudo que é multiplicado por zero da zero, mas temos uma reta, só que sem inclinação, ou melhor, com inclinação zero) e depois do 5 temos uma reta f(x)=1+2x.

Nesse caso, o partição, o valor 5 tem um significa biológico importante, ele é o limite mínimo de nutrientes, para começarmos a observar alguma coisa, antes dele nada acontece, mas como podemos recuperar esse valor dos dados? A primeira opção é usar um tipo de força bruta, testar o modelo acima.

Para aquela equação podemos ajustar ela usando uma variável tipo dummy usando a função lm, ou usar a função nls e ajustar diretamente a equação piecewise usando um ifelse, aqui vamos usar diretamente a equação com nls que é mais direto a interpretação dos parâmetros(uma reta é a+bx, então vamos ter dois a, o a1 e a2 e dois b, o b1 e b2 para as duas retas), ajustamos duas retas e um valor de quebra, no caso ajustamos o modelo para muitas valores de quebra e vemos qual da o melhor ajuste usando logLikelihood

Então criamos um vetor de possíveis valores de quebra.

1
2
quebras <- seq(1,9,0.5)
minimo_quadrado <- vector(length = length(quebras))

E testamos todos eles

1
2
3
4
for(i in 1:length(quebras)){
 piecewise <- nls(tamanho ~ ifelse(nutriente < quebras[i],a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1))
 minimo_quadrado[i] <- logLik(piecewise)
}

E podemos avaliar graficamente o resultado.

Certo, como esperávamos, o valor da simulação é recuperado, veja que o pico de ajusta está bem no 5.

Podemos então ajustar o modelo usando 5 como ponto de quebra e ver se recuperamos os coeficientes corretamente

> modelo<-nls(tamanho ~ ifelse(nutriente < 5,a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1)) > summary(modelo) Formula: tamanho ~ ifelse(nutriente < 5, a1 + b1 * nutriente, a2 + b2 * nutriente) Parameters: Estimate Std. Error t value Pr(>|t|) a1 7.7307 0.3410 22.668 <2e-16 *** a2 3.2449 1.1431 2.839 0.0063 ** b1 0.2193 0.1249 1.756 0.0846 . b2 1.7430 0.1369 12.734 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9221 on 56 degrees of freedom Number of iterations to convergence: 1 Achieved convergence tolerance: 1.57e-08

E temos as inclinações, zero para a primeira equação esta no intervalo do b1, que é o nome que demos para o parâmetro, e assim como o 2 está no intervalo do b2, os interceptos também são recuperados. Vamos ver ele graficamente, e além disso, vamos comparar com um modelo linear.

Veja que ele fica bem mais certinho, mas o modelo linear também detectaria uma tendência aqui, porém não nos forneceria essa informação, de a partir de onde, de qual concentração é relevante, seriamos enganados pensando que a concentração de nutrientes é relevante em qualquer intensidade quando ela so importa em uma parte do preditor, isso mesmo avaliando os resíduos dos modelos.

Temos um ajuste menor (lembre-se que estamos falando de números negativos, então -78 é maior que -111), mas não temos uma tendência clara nos resíduos para esse caso, rejeitar o modelo linear, somente observando os resíduos, então aqui interessa muito o que pensamos a priori da construção do modelo.

Como podemos ter casos onde a busca precisa ser mais precisa para o ponto de quebra, temos funções prontas que vão cuidar dessas atividades, por exemplo o pacote segmented, que ajusta esse tipo de modelo diretamente, fazendo um update do modelo linear.

Basicamente temos o seguinte resultado:

> modelo_pieciwise<- segmented(modelo_linear, seg.Z = ~nutriente, psi=1) > summary(modelo_pieciwise) ***Regression Model with Segmented Relationship(s)*** Call: segmented.lm(obj = modelo_linear, seg.Z = ~nutriente, psi = 1) Estimated Break-Point(s): Est. St.Err 3.799 0.304 Meaningful coefficients of the linear terms: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.8733 0.4533 17.369 <2e-16 *** nutriente 0.1274 0.1967 0.648 0.52 U1.nutriente 1.9158 0.2222 8.621 NA --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.045 on 56 degrees of freedom Multiple R-Squared: 0.9576, Adjusted R-squared: 0.9554 Convergence attained in 7 iterations with relative change 3.485019e-16

Veja que temos o Estimated Break-Point, que da 3.799, um pouco inferior ao que realmente usamos, e temos apenas um intercepto, mas a primeira inclinação fica em zero, e a segunda em 0.12+1.91, diferença de inclinação para depois do ponto de quebra, que é 2, então ele estima bem os parâmetros, mas é um modelo um pouco diferente do que usamos (só um intercepto), eu não li muito sobre esse pacote, para saber ajustar exatamente como fizemos ali em cima, mas é um pacote a ser explorado, para quem quer utilizar esse tipo de modelo e o sistema de plot dele é bem simples.

Bem é isso ai, depois de um tempo enroscado, de volta a atividade, o script vai estar lá no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Michael J. Crawley 2012 The R Book, 2nd ed. Wiley 1076pp

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
set.seed(31415)
 
nutriente <- runif(60,0,10)
 
fn_bs <- function(x){
    saida<-vector(length=length(x))
    for(i in 1:length(x)){       
        if(x[i]<5){
            saida[i] <- 8+0*x[i]
        }else{
            saida[i] <- 1+2*x[i]
        }
    }
    saida <- saida+rnorm(length(x),0,1)
    return(saida)
}
 
tamanho <- fn_bs(nutriente)
 
plot(tamanho~nutriente,xlab="Concentração de nutrientes",ylab="Tamanho da plantula (cm)",frame=F,ylim=c(0,30),xlim=c(0,12),pch=19)
 
#####################################
## Força bruta
#####################################
quebras <- seq(1,9,0.5)
minimo_quadrado <- vector(length = length(quebras))
 
for(i in 1:length(quebras)){
 piecewise <- nls(tamanho ~ ifelse(nutriente < quebras[i],a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1))
 minimo_quadrado[i] <- logLik(piecewise)
}
 
minimo_quadrado <- unlist(minimo_quadrado)
 
plot(quebras,minimo_quadrado,type="b",pch=19,xlab="Valor de quebra",ylab="LogLikelihood",frame=F,xlim=c(0,10))
 
modelo<-nls(tamanho ~ ifelse(nutriente < 5,a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1))
summary(modelo)
 
modelo_linear<- lm(tamanho~nutriente)
summary(modelo_linear)
 
plot(tamanho~nutriente,xlab="Concentração de nutrientes",ylab="Tamanho da plantula (cm)",frame=F,ylim=c(0,30),xlim=c(0,12),pch=19)
 
abline(modelo_linear,lwd=3,lty=3,col="blue")
 
curve(coef(modelo)[1]+coef(modelo)[3]*x,0,5,add=T,lwd=3,lty=3,col="red")
curve(coef(modelo)[2]+coef(modelo)[4]*x,5,12,add=T,lwd=3,lty=3,col="red")
abline(v=5,lwd=1,lty=2,col="red")
legend("topleft",lwd=3,lty=3,col=c("blue","red"),legend=c("Modelo Linear","Modelo piecewise"),bty="n")
 
par(mfrow=c(2,1))
plot(resid(modelo_linear),main=paste("Loglikelihood do modelo liner =",round(logLik(modelo_linear),3)),frame=F,xlab="Amostras",ylab="Resíduo",ylim=c(-5,5))
abline(h=0,lty=2,lwd=2)
plot(resid(modelo),main=paste("Loglikelihood do modelo piecewise =",round(logLik(modelo),3)),frame=F,xlab="Amostras",ylab="Resíduo",ylim=c(-5,5))
abline(h=0,lty=2,lwd=2)
 
#####################################
## Pacote Segmented
#####################################
##install.packages("segmented")
library(segmented)
 
modelo_pieciwise<- segmented(modelo_linear, seg.Z = ~nutriente, psi=1)
summary(modelo_pieciwise)
 
plot(tamanho~nutriente,xlab="Concentração de nutrientes",ylab="Tamanho da plantula (cm)",frame=F,ylim=c(0,30),xlim=c(0,12),pch=19)
plot(modelo_pieciwise,add=T,lwd=2,lty=2,col="red")