Writing a Monte Carlo simulation in R I am trying to write a Monte Carlo simulation in R and I am really stuck! I want to know the probability distribution of a random person in the UK becoming ill from eating a cooked 100g piece of chicken. 
I have the following information: out of 1000 pieces of chicken tested 20 had bacteria in question and I have data for the $\log_{10}$ counts of these 20 pieces, I also have min and max $\log_{10}$ counts (0.1 and 3.0). I also know the average person in UK eats 2 x 100g portions of chicken a week. The model for risk of illness given an ingested number of the bacteria is predicted by $R=1-\exp(-aD)$ where $D$ is the ingested number of organisms and I have a value for $a$.
I can write basic Monte Carlo simulations but I am struggling with the start of this one as I can't get my head around the model being ingested bacteria and the question being risk from eating a 100g portion.


*

*Is my first step here to obtain the CDF? 

*And what is the distribution I should use?

 A: Here my 2 cents.  When 20/1,000 samples are contaminated, you could use a Poisson distribution to model the event of a person eating a contaminated sample (mind the difference between the weight of a sample in the lab and on your plate).
When a Monte Carlo sample is simulated to be positive in the previous step, you can then simulate the ingested number of bacteria using the empirical distribution of your 20 values, or fit for example an exponential distribution.  In the latter case, have a look at the fitdistrplus package.  In the former case, you could do it as follows:
Fn <- ecdf(d)
random_samples <- quantile(x=Fn, prob=runif(n=1e05, min=0, max=1))

Take care that you transfer your units correctly (log values, sample weights, etc.) and realize that this simulation is only a very crude approximation of what happens in reality.
A: Here my 2 cents as well! With the given information, the probability distribution of D can be derived as follows:
Pr(D)=Pr(D|Contaminated)*P(Contaminated)

We can estimate Pr(Contaminated)=20/1000=1/50. For Pr(D|Contaminated) a simple histogram suggests a Poisson like distribution. But you can fit another distribution. However, for count data it is customary to use Poisson or Binomial distributions. 
Once you have the distribution for D, you can sample from it and plugin in the risk of illness, R(D). This risk can be interpreted as the probability of been sick in one trial (i.e. eating a piece). Then the distribution of been sick follows S ~ Binomial(R(D),2) if each piece is 100g. The probability for one person to be sick in one week is then Pr(S=1)=2*R(D)*[1-R(D)].
To summarize:


*

*Fit a distribution to your D s (10^d from you comments). Possibly a Poisson.

*Sample a value from this distribution (truncate if necessary to take into account the minimum and the maximum).

*Compute 2*R(D)*[1-R(D)]

*Go to step 2

