def nb_fit(data, P_init=None, R_init=None, epsilon=<NUM_LIT>, max_iters=<NUM_LIT:100>):
means = data.mean(<NUM_LIT:1>)<EOL>variances = data.var(<NUM_LIT:1>)<EOL>if (means > variances).any():<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>genes, cells = data.shape<EOL>P = <NUM_LIT:1.0> - means/variances<EOL>R = means*(<NUM_LIT:1>-P)/P<EOL>for i in range(genes):<EOL><INDENT>result = minimize(nb_ll_row, [P[i], R[i]], args=(data[i,:],),<EOL>bounds = [(<NUM_LIT:0>, <NUM_LIT:1>), (eps, None)])<EOL>params = result.x<EOL>P[i] = params[<NUM_LIT:0>]<EOL>R[i] = params[<NUM_LIT:1>]<EOL><DEDENT>return P,R<EOL>
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