Non-Poissonian satellite occupation model

[aphearin]

This PR brings in a model for non-Poissonian occupation statistics of satellites. The new satellite occupation model MBK10Sats implements a clever technique due to @mrbk to model satellite HODs based on the Negative Binomial distribution. Deviations from Poisson fluctuations of satellites should generically be expected at some level since it is commonplace for subhalos to accrete in groups (aka "satellite preprocessing"), and recent papers have explored this effect within the context of HOD modeling using somewhat ad hoc methods.

In the No Place Like Home paper, the Negative Binomial distribution was shown to provide an accurate characterization of the second occupation moment of subhalos of Milky Way mass hosts in the MS-II and Aquarius simulations. This distribution has a nice implementation in scipy that enables both analytical calculations and Monte Carlo realizations to be straightforward and fast.

The MBK10Sats class has an identical first moment as Zheng07Sats (the usual power law function of halo mass). Whereas the variance in a Poisson system is just equal to the first moment, the second moment in the MBK10Sats model has an additional dependence upon a single new parameter, nsat_up0, which controls the deviation from Poisson fluctuations. The MBK10Sats model behaves as a Poisson model in the limit of infinite nsat_up0; smaller values of nsat_up0 produce stronger deviations from Poisson fluctuations (with larger variance relative to Poisson).

Specifically, the variance in a Negative Binomial distribution distribution is controlled by a parameter 0 < p < 1, and in the implementation brought in by this PR, the parameter nsat_up0 is related to p via a simple sigmoid transformation. This way, in an MCMC, when the walkers explore the parameter space, the nsat_up0 parameter can take on any value without raising an exception that halts the burn. The non_poissonian_p method of MBK10Sats class can be used to provide the value of p corresponding to the nsat_up0 parameter. The figure below shows how the strength of Poisson deviations depend upon the nsat_up0 parameter.

This next figure reproduces Figure 10 of No Place Like Home by using the MBK10Sats class to populate a sample of massive halos.

This PR also brings in the AssembiasMBK10Sats class, in which the mean number of satellites exhibits assembly bias at fixed mass, where the assembly bias is captured with the Decorated HOD. The AssembiasMBK10Sats class inherits from the AssembiasZheng07Sats class, and overloads the second moment using the Negative Binomial distribution.

The convention adoped in the AssembiasMBK10Sats class is that P(Nsat | M, x) is modeled using the Negative Binomial distribution, where x is the secondary halo property that controls galaxy assembly bias such that <Nsat | M, x> != <Nsat | M>. This is in contrast to AssembiasZheng07Sats, which models P(Nsat | M, x) according to a Poisson distribution. Note that in both models, the occupation statistics P(Nsat | M) is non-Poissonian, which is a generic feature of any assembly bias model that preserves the first moment (see section 4.1.3 of the Decorated HOD paper for further discussion of the relationship between assembly bias and the second satellite moment).