Rosalind – Enumerating k-mers Lexicographically

Esse problema é muito simples, dada um alfabeto, uma coleção com um certo número de símbolos e um tamanho, temos que fazer todas as combinações possíveis organizados lexicograficamente, na “ordem alfabética”.

Primeiro, a quantidade de possibilidades será o número de símbolos elevado pelo tamanho da palavra que vamos formar.

O exemplo de dados é:

A C G T 2

Então o alfabeto é o ‘A’, ‘C’, ‘G’ e ‘T’ e o tamanho é 2, ou seja temos 2 posições

_ _

E cada posição tem quatro possibilidades.

A primeira solução que eu pensaria é fazer um loop, usando um for para cada posição, algo como:

1
2
3
4
5
alfabeto=['a','c','g','t']
 
for i in alfabeto:
    for j in alfabeto:
        print i+j

Mas ai temos o problema de como criaríamos for de acordo com o tamanho de entrada, já que ele é variável. Mas podemos criar uma solução recursiva para esse problema.

1
2
3
4
5
6
7
8
def fn_re(n):
 
    if n>1:        
        fn_re(n-1)
 
    print n
 
fn_re(len(alfabeto))
>>>4 3 2 1

Ou seja, enquanto n for maior que o tamanho for maior que zero, a gente faz algo, até a condição de parada, que é n chegar a zero, não ter o que fazer, ai a gente só retorna o processamento. Então o que a gente precisa é mandar o alfabeto, e mandar o tamanho da palavra que vamos formar, ai para cada letra do alfabeto, a gente chama de novo a função.

1
2
3
4
5
6
7
def combina(alfabeto,n,kmer,lista):
    if n==0:
        lista.append(kmer)
    else:
        for letra in alfabeto:
            combina(alfabeto,n-1,kmer+letra,lista)
    return lista

Um jeito legal de entender melhor o que está acontecendo, é imprimir o kmer e a lista, no inicio da função, assim:

1
2
3
4
5
6
7
8
9
def combina(alfabeto,n,kmer,lista):
    print kmer
    print lista
    if n==0:
        lista.append(kmer)
    else:
        for letra in alfabeto:
            combina(alfabeto,n-1,kmer+letra,lista)
    return lista

Ai quando a gente usa ela

1
2
3
alfabeto=['a','c','g','t']
n=2
resposta= combina(alfabeto,n,'',[])

A gente vê o seguinte:

[] a [] aa [] ac ['aa'] ag ['aa', 'ac'] at ['aa', 'ac', 'ag'] c ['aa', 'ac', 'ag', 'at'] ca ['aa', 'ac', 'ag', 'at'] cc ['aa', 'ac', 'ag', 'at', 'ca'] . . .

A primeira chamada, o kmer é ” e a lista ta vazia.
O n é diferente de zero, então a gente continua e então chamamos a função para todo o alfabeto, no primeira chamada recursiva, vamos diminuir um no n, que vai ficar 1 e veja que vamos concatenar no for a uma letra ao kmer, a primeira letra do nosso alfabeto é o a, então na segunda vez que imprimimos na tela, vemos o a, porque fizemos kmer+letra na chamada e temos o n sendo 1, então chamamos recursivo de novo, só que agora, na recursão o kmer é ‘a’, e no for vamos fazer kmer+letra de novo, e chamar recursivo, mas agora temos ‘aa’, agora a gente vai cair dentro do if, e ao invés de fazer a chamada recursiva, a gente so coloca o ‘aa’ na lista, e assim seguimos.

Bem agora é só ler o arquivo, pegar os parâmetros e pimba, 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.

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
def combina(alfabeto,n,kmer,lista):
    print kmer
    print lista
    if n==0:
        lista.append(kmer)
    else:
        for letra in alfabeto:
            combina(alfabeto,n-1,kmer+letra,lista)
    return lista
 
"""
file = open( "/mnt/B21AA1BD1AA17ECB/downloads_chromium_ubuntu/rosalind_lexf.txt" )
alfabeto = file.readline().split()
n = int(file.readline())
file.close()
"""
 
###Exemplo
alfabeto=['a','c','g','t']
n=2
 
 
print "Número de possibilidades: "+str(len(alfabeto)**n)
 
resposta= combina(alfabeto,n,'',[])
 
for i in resposta:
    print i

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)