How to sample for conditional probability from unknown populations 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.
 A: 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

A: 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)
