The parameter is the unknown species, either A or B. Prior probabilities are equal (to 1/2). The datum is the observation of a twin birth, and the likelihood function is the probability of a twin birth for each value of p. The posterior distribution is the normalized product of the two:
> prior1 <- rep(0.5, 2) # Pr(p = A) = Pr(p = B) = 1/2
> names(prior1) <- c("A", "B")
> likelihood1 <- c(0.1, 0.2) # Pr(twin | p) for p = A, B
> names(likelihood1) <- c("A", "B")
> posterior1 <- prior1*likelihood1
> library(MASS) # for fractions()
> fractions(posterior1 <- posterior1/sum(posterior1)) # Pr(p | twin), normalized
A B
1/3 2/3
Then by the law of total probability:
> fractions(sum(posterior1*likelihood1)) # Pr(p = A) Pr(twin | p = A) + Pr(p = B) Pr(twin | p = B)
[1] 1/6
The posterior probability of species A was already computed:
> fractions(posterior1["A"])
A
1/3
The updated prior is the previous posterior, and the new datum is the observation of a singleton birth:
> prior2 <- posterior1
> likelihood2 <- 1 - likelihood1 # Pr(singleton | p) = 1 - Pr(twin | p)
> posterior2 <- prior2*likelihood2
> (posterior2 <- posterior2/sum(posterior2))["A"]
A
0.36
Ignoring the information on births returns us to the equal-probability prior for species A or B. The datum is the test result for A.
> likelihood3 <- c(0.8, 1 - 0.65) # Pr (test is A | p) for p = A, B
> names(likelihood3) <- c("A", "B")
> posterior3 <- prior1*likelihood3
> (posterior3 <- posterior3/sum(posterior3))["A"]
A
0.6956522
Using the information on births, we take the posterior after the second birth (i.e., twin followed by singleton) as the prior:
> prior3 <- posterior2
> posterior4 <- prior3*likelihood3
> (posterior4 <- posterior4/sum(posterior4))["A"]
A
0.5625