How to use the lifetimes.generate_data.beta_geometric_beta_binom_model function in Lifetimes

def test_params_same_from_sim_data(self, bbgb_params):
        sim_data = beta_geometric_beta_binom_model(N=6, size=100000, **bbgb_params)
        bbtf = estimation.BetaGeoBetaBinomFitter()
        grouped_data = sim_data.groupby(['frequency', 'recency', 'n_periods'])['customer_id'].count()
        grouped_data = grouped_data.reset_index().rename(columns={'customer_id': 'weights'})
        bbtf.fit(grouped_data['frequency'],
                 grouped_data['recency'],
                 grouped_data['n_periods'],
                 grouped_data['weights'])

        npt.assert_allclose(
            np.asarray(list(bbgb_params.values())).astype(float),
            np.asarray(bbtf._unload_params('alpha', 'beta', 'gamma', 'delta')).astype(float),
            atol=0.1, rtol=1e-2)

def test_hitting_max(self, bbgb_params):
        sim_data = beta_geometric_beta_binom_model(N=6, size=5000, **bbgb_params)
        assert sim_data['frequency'].max() == 6
        assert sim_data['recency'].max() == 6

def test_positivity(self, bbgb_params):
        sim_data = beta_geometric_beta_binom_model(N=6, size=5000, **bbgb_params)
        assert (sim_data['frequency'] >= 0).all()
        assert (sim_data['recency'] >= 0).all()

def test_alive_probs(self, bbgb_params):
        sim_data = beta_geometric_beta_binom_model(N=6, size=50000, **bbgb_params)
        assert (np.abs(sim_data.loc[(sim_data['frequency'] == 0) & (sim_data['recency'] == 0),
                                    'alive'].mean() - 0.11) < 0.01)
        assert (np.abs(sim_data.loc[(sim_data['frequency'] == 2) & (sim_data['recency'] == 4),
                                    'alive'].mean() - 0.59) < 0.01)
        assert (np.abs(sim_data.loc[(sim_data['frequency'] == 6) & (sim_data['recency'] == 6),
                                    'alive'].mean() - 0.93) < 0.01)

        self.generate_new_data = lambda size=1: beta_geometric_beta_binom_model(
            # Making a large array replicating n by n_custs having n.
            np.array(sum([n_] * n_cust for (n_, n_cust) in zip(n_periods, weights))),
            *self._unload_params("alpha", "beta", "gamma", "delta"),
            size=size
        )