I'm trying to generate "species counts" in raster form -- essentially, the number of individuals in a species in a raster cell. I'd like to use the negative binomial distribution to do so. Normally I would just use something like the Create Random Raster tool in ArcGIS; however, there is an extra wrinkle. I'd also like the species to be weighted based on a second raster, where higher values are more likely to fall in certain areas.
So if I were hypothetically modeling a species like a lion, I'd create a raster with possible values taken from a negative binomial distribution where higher values were more likely on savannahs, less likely in urban areas, and never found in the middle of lakes or oceans (with the land classes taken from an NDVI raster, or something similar).
I'm not sure how to do this, and would really appreciate any pointers on ArcGIS tools, R packages/functions, etc.