1
$\begingroup$

I am providing the full question as well my solution below. I'm looking for help with part (d), a simulation question.

Q - Suppose there are two species of Pandas, $T_1$ and $T_2$ which are indistinguishable and exist in equal proportions, but differ in how they lay children. Species $T_1$ gives birth to twins 10% of the time and otherwise lays a single cub. Species $T_2$ lays twins 20% of the time, and otherwise only lays a single cub.

There are two pandas, who are unrelated and of unknown species, Panda X and Panda Y.

a) Panda X has twins the first year. Find the conditional probability that the Panda is species $T_1$.

b) Find the conditional probability that Panda X will have twins again.

c) Suppose there exists a genetic test which correctly identifies pandas as species $T_1$ 80% of the time and correctly identifies species $T_2$ 60% of the time. This test is administered to Panda Y, and the results indicate that the panda is species $T_1$. Find the probability that the first birth from Panda Y results in twins.

d) Verify this through simulation in R.

Answers: (a) and (b) I do not have issues with. Providing part (c) for reference.

c) Follows Bayes laws again, but it is a new panda so we can forget parts a) and b). We use the Bayes rule to get the posterior probability they are type $S_1$ and then use the approach from part (b) to compute the probability the birth is twins. For the first part,

$$P(S_1\mid test) = \frac{P(test\mid S_1)\frac{1}{2}}{P(test\mid S_1)\frac{1}{2} + P(test\mid S_2)\frac{1}{2}} = \frac{\frac{8}{10}\frac{1}{2}}{\frac{8}{10}\frac{1}{2}+\frac{4}{10}\frac{1}{2}} = \frac{{\frac{2}{5}}}{\frac{3}{5}} = \frac{2}{3}$$

since the probability that it is right given it's a $T_1$ is $8/10$ and the probability the test is wrong given that it's a $T_2$ is $4/10.$

My question is, how do I go about part d? I'm aware that I'd have to sample in some way, but considering I don't know how the populations are distributed, I can't take any samples from pre-existing functions within R. Would be great if someone could show me how to solve part d. Many thanks in advance.

Improve this question
$\endgroup$
3
  • $\begingroup$ Answer (c) is incorrect. And to run a simulation confirmation of the results, you have to generate a large number of pandas, by first generating species (binary), presence of twins (binary), second birth presence of twins (binary), detection of species T1 (binary). [Actually, imho, I find the wording of the exercise rather poor, including the lack of the event against which to condition in (b).] $\endgroup$
    – Xi'an
    Commented Oct 25, 2018 at 4:44
  • $\begingroup$ Is the methodology for (c) incorrect or simply the answer itself? I get what you're saying about the generating sample part, but it's the coding itself I'm struggling with... $\endgroup$
    – user223357
    Commented Oct 25, 2018 at 4:50
  • $\begingroup$ The probability of having twins is necessarily between 10% and 20%. Generating a binary variable in R can be done by sample or alternatively by generating a uniform by runif and comparing it to the probability. $\endgroup$
    – Xi'an
    Commented Oct 25, 2018 at 5:01

2 Answers 2

Reset to default
0
$\begingroup$

Xi'an's answer gives good general pseudo-code for such a problem, but I thought it might also be helpful to have someone tackle the R-specific aspect of this, which your post kind of implies you were having some trouble with. You seemed skeptical that you could use the builtin sample functions from R, but you surely can:

N <- 1e5 # Total number of pandas in simulation -- trying with 100K
set.seed(42) # Set the seed for reproducibility
types <- sample(x = 1:2, size = N, replace = TRUE) # Generate types
twins <- sapply(types, function(x) { # Was the first litter twins?
    sample(x = 0:1, size = 1, prob = c(1 - x*0.1, x*0.1))
})
test_results <- sapply(types, function(x) { # What did the type test say?
    p <- 0.8*(2-x) + 0.4*(x-1)
    sample(x = 1:2, size = 1, prob = c(p, 1 - p))
})
tested_as_T1 <- which(test_results == 1) # Which ones tested as type 1?
# And the finale.... What proportion of those had twins?
sum(twins[tested_as_T1]) / length(tested_as_T1)
#> [1] 0.130611
# What was the theoretical answer again?
(2/3) * 0.1 + (1/3) * 0.2
#> [1] 0.1333333

Created on 2018-10-25 by the reprex package (v0.2.1)

Improve this answer
$\endgroup$
2
  • $\begingroup$ hmmm.... I refrained from providing the complete R code since this was a self-study question where the OP is asking for help in solving an exercise rather than having someone else doing it. $\endgroup$
    – Xi'an
    Commented Oct 26, 2018 at 8:41
  • $\begingroup$ @Xi'an I see. I think we differently interpreted which part(s) of the question OP wanted more hand-holding on. I had thought the "self study" part was the stats aspect, i.e. that they are trying to figure out the stats parts on their own, but wanted more explicit help in how to get a simulation going in R, but in retrospect I can totally see your interpretation as well. $\endgroup$
    – duckmayr
    Commented Oct 26, 2018 at 8:52
1
$\begingroup$

A Monte Carlo code to check the first three answers would look like

for p in {1,...,P} do 
  simulate and store type t[p]
  simulate and store first birth size f[p]
  simulate and store second birth size s[p]
  simulate and store test result r[p]
endfor

countt=0
counts=0
countf
for p such that f[p]=2
  countf=countf+1
  countt=countt+1 if t[p]=1
  counts=counts+1 if s[p]=2
endfor
countt=countt/countf
counts=counts/countf

countr=0
countt=0
for p such that r[p]=1
   countr=countr+1
   countt=countt+1 if f[p]=1
endfor
countt=countt/countr
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.