We vinden op straat een munt en willen graag weten of deze zuiver is.

We gooien 10 keer, waarvan 2 keer kop. De binomiaal verdeling hiervan ziet er als volgt uit.

k = 10
n = 0:k # Outcome spece

head = 2

d1 = dbinom(n, k, head/k)


title = "Binomiaal verdelingen"

colh0 = "orange"
colha = "green"

barplot(d1,
        main        = title, 
        names.arg   = 0:k, 
        xlab        = "Aantal kop", 
        ylab        = "Percentage", 
        ylim        = c(0,max(d1)+.1),
        col         = rbind(colha),
        #col        = "white",        
        beside      = TRUE,
        legend.text = c("H1"),
        args.legend = list(x = "topleft")
        ) -> bp

alt = c(rbind(d1))
text(bp, alt, round(alt, 2),
     #pos = 3,
     adj = c(-.3,.5),
     cex = .75, 
     srt = 90)

Wat is de “likelyhood” voor elke mogelijke uitkomst (0 tot en met 10), afgezet tegen de \(H_0\) verwachting dat we 5 keer kop zullen gooien?

Om hier achter te komen moeten we kijken naar de verhouding van de verwachte uitkomst onder \(H_1\) \(P(D|H_1)\) met de verwachte kans op 5 keer kop onder \(H_1\).

# Probability data given h1
pdh1 <- round(alt, 2)
pdh1
##  [1] 0.11 0.27 0.30 0.20 0.09 0.03 0.01 0.00 0.00 0.00 0.00
lh <- pdh1 / 0.03
lh
##  [1]  3.6666667  9.0000000 10.0000000  6.6666667  3.0000000  1.0000000
##  [7]  0.3333333  0.0000000  0.0000000  0.0000000  0.0000000

Om te komen tot een gewogen “likelyhood” (WLH) moeten we de likelyhood vermenigvuldigen met onze prior.

Laten we als prior kiezen voor een \(H_0\) binomiaal verdeling met een verwachte waarde van .5.

h0 = .5

d0 = dbinom(n, k, h0)
barplot(rbind(d1, d0),
        main        = title, 
        names.arg   = 0:k, 
        xlab        = "Aantal kop", 
        ylab        = "Percentage", 
        ylim        = c(0,max(c(d1, d0))+.1),
        col         = rbind(colha,colh0),
        #col        = "white",        
        beside      = TRUE,
        legend.text = c("H1", "Prior"),
        args.legend = list(x = "topleft")
        ) -> bp

alt = c(rbind(d1, d0))
text(bp, alt, round(alt, 2),
     #pos = 3,
     adj = c(-.3,.5),
     cex = .75, 
     srt = 90)

De WLH komt tot stand door de lh te vermenigvuldigen met de kansen op elke afzonderlijke uitkomst onder onze prior, in dit geval \(H_0\).

lh
##  [1]  3.6666667  9.0000000 10.0000000  6.6666667  3.0000000  1.0000000
##  [7]  0.3333333  0.0000000  0.0000000  0.0000000  0.0000000
d0
##  [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 0.2050781250
##  [6] 0.2460937500 0.2050781250 0.1171875000 0.0439453125 0.0097656250
## [11] 0.0009765625
wlh = lh * d0
wlh
##  [1] 0.003580729 0.087890625 0.439453125 0.781250000 0.615234375
##  [6] 0.246093750 0.068359375 0.000000000 0.000000000 0.000000000
## [11] 0.000000000

Wat ons dat oplevert is de posterior verdeling.

barplot(wlh,
        main        = "Posterior", 
        names.arg   = 0:k, 
        xlab        = "Aantal kop", 
        ylab        = "Percentage", 
        ylim        = c(0,max(wlh)+.1),
        col         = "blue",
        #col        = "white",        
        beside      = TRUE,
        legend.text = c("Posterior"),
        args.legend = list(x = "topleft")
        ) -> bp

alt = c(rbind(d1))
text(bp, wlh, round(wlh, 2),
     #pos = 3,
     adj = c(-.3,.5),
     cex = .75, 
     srt = 90)

De Bayes factor voor \(H_1\) is het gemiddelde van de posterior verdeling. Met andere woorden het gemiddelde van de gewogen likelyhood (wlh).

\(\text{BF}_{10}\)

mean(wlh)
## [1] 0.2038056

Wat duidelijk moet zijn is dat je dus geen alternatieve hypothese specificeerd maar dat je kijkt wat de verdeling is gegeven de data en dat afzet tegen je aanvankelijke verwachting. In dit geval een \(H_0\) verdeling.

Volgens mij klopt onderstaande niet

In termen van Bayes theorem

\[P(A\mid B) = \frac{P(B \mid A) \, P(A)}{P(B)}\]

vermenigvuldigen we de likelyhood \(P(B \mid A)\) met de prior \(P(A)\) om te komen tot onze posterior verdeling. De kans op de data \(P(B)\) kunnen we in deze achterwege laten aangezien dit weer een constante is in de Bayes factor.

Nog een keer

Stel nu dat we nog 10 keer gooien met de munt en we dit maal 6 keer kop gooien. Wat voor invloed heeft dat dan op de posterior verdeling en de Bayes factor?

Onze nieuwe prior wordt de vorige posterior

h0 = 2/10
d0 = wlh


d1 = dbinom(n,k,6/10)

barplot(d1,
        main        = title, 
        names.arg   = 0:k, 
        xlab        = "Aantal kop", 
        ylab        = "Percentage", 
        ylim        = c(0,max(c(d1, d0))+.1),
        col         = rbind(colha),
        #col        = "white",        
        beside      = TRUE,
        legend.text = c("H1"),
        args.legend = list(x = "topleft")
        ) -> bp

alt = c(rbind(d1))
text(bp, alt, round(alt, 2),
     #pos = 3,
     adj = c(-.3,.5),
     cex = .75, 
     srt = 90)

# Probability data given h1
pdh1 <- round(alt, 2)
pdh1
##  [1] 0.00 0.00 0.01 0.04 0.11 0.20 0.25 0.21 0.12 0.04 0.01
lh <- pdh1 / 0.01
lh
##  [1]  0  0  1  4 11 20 25 21 12  4  1
barplot(rbind(d1, d0),
        main        = title, 
        names.arg   = 0:k, 
        xlab        = "Aantal kop", 
        ylab        = "Percentage", 
        ylim        = c(0,max(c(d1, d0))+.1),
        col         = rbind(colha,colh0),
        #col        = "white",        
        beside      = TRUE,
        legend.text = c("H1", "Prior"),
        args.legend = list(x = "topleft")
        ) -> bp

alt = c(rbind(d1, d0))
text(bp, alt, round(alt, 2),
     #pos = 3,
     adj = c(-.3,.5),
     cex = .75, 
     srt = 90)

lh
##  [1]  0  0  1  4 11 20 25 21 12  4  1
d0
##  [1] 0.003580729 0.087890625 0.439453125 0.781250000 0.615234375
##  [6] 0.246093750 0.068359375 0.000000000 0.000000000 0.000000000
## [11] 0.000000000
wlh = lh * d0
wlh
##  [1] 0.0000000 0.0000000 0.4394531 3.1250000 6.7675781 4.9218750 1.7089844
##  [8] 0.0000000 0.0000000 0.0000000 0.0000000
barplot(wlh,
        main        = "Posterior", 
        names.arg   = 0:k, 
        xlab        = "Aantal kop", 
        ylab        = "Percentage", 
        ylim        = c(0,max(wlh)+.1),
        col         = "blue",
        #col        = "white",        
        beside      = TRUE,
        legend.text = c("Posterior"),
        args.legend = list(x = "topleft")
        ) -> bp

alt = c(rbind(d1))
text(bp, wlh, round(wlh, 2),
     #pos = 3,
     adj = c(-.3,.5),
     cex = .75, 
     srt = 90)