def fit_mme(sample_data, replace_values=None, name=None):
"""
    Calculates the Maximum-of-Momentum Estimator of `NB(r, p)` for given sample data along axis 0.

    :param sample_data: matrix containing samples for each distribution on axis 0n
        E.g. `(N, M)` matrix with `M` distributions containing `N` observed values each
    :param replace_values: Matrix of size `shape(sample_data)[1:]`
    :param name: A name for the operation (optional).
    :return: estimated values of `r` and `p`
    """
    with tf.name_scope(name, "MME"):
        mean = tf.reduce_mean(sample_data, axis=0, name="mean")
        variance = tf.reduce_mean(tf.square(sample_data - mean),
                              axis=0,
                              name="variance")
        if replace_values is None:
            replace_values = tf.fill(tf.shape(variance), math.nan, name="NaN_constant")

        r_by_mean = tf.where(tf.less(mean, variance),
                         mean / (variance - mean),
                         replace_values)
        r = r_by_mean * mean
        r = tf.identity(r, "r")

        p = 1 / (r_by_mean + 1)
        p = tf.identity(p, "p")

        return r, p