I have count data that is bounded on one side at zero (see image). It is bimodal and I think it results from two different processes. I would like to fit a poisson distribution around the hump around zero, and a gaussian distribution around the hump at higher values. Then, I would like to find a cutoff separating the values into the two distributions. However, I'm struggling to find a method in R to do so. Does anyone know of an R package that can do this?