<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def nb_fit(data, P_init=None, R_init=None, epsilon=1e-8, max_iters=100):
 """ Fits the NB distribution to data using method of moments. Args: data (array):
 genes x cells P_init (array, optional):
 NB success prob param - genes x 1 R_init (array, optional):
 NB stopping param - genes x 1 Returns: P, R - fit to data """

    means = data.mean(1)
    variances = data.var(1)
    if (means > variances).any():
        raise ValueError("For NB fit, means must be less than variances")
    genes, cells = data.shape
    # method of moments
    P = 1.0 - means/variances
    R = means*(1-P)/P
    for i in range(genes):
        result = minimize(nb_ll_row, [P[i], R[i]], args=(data[i,:],),
                bounds = [(0, 1), (eps, None)])
        params = result.x
        P[i] = params[0]
        R[i] = params[1]
        #R[i] = fsolve(nb_r_deriv, R[i], args = (data[i,:],))
        #P[i] = data[i,:].mean()/(data[i,:].mean() + R[i])
    return P,R