A sample of bay anchovy larvae counts yielded the data

158, 143, 106, 57, 97, 80, 109, 109, 350, 224, 109, 214, 84

The negative binomial(r,p) distribution (with mean m= r(1-p)/p and
variance m + (m^2)/r ) is a reasonable model. For m=xbar, what values of r lead to simulated distributions that are consistent with the data above?

Simulation should be done with SAS or R.