xs = np.random.uniform(low, high, size=100)
means = intercept + slope * xs
ys = np.array([stats.poisson(mu).rvs() for mu in means])