def estimate_row_means(
            self,
            X,
            observed,
            column_means,
            column_scales):
        """
        row_center[i] =
        sum{j in observed[i, :]}{
            (1 / column_scale[j]) * (X[i, j] - column_center[j])
        }
        ------------------------------------------------------------
        sum{j in observed[i, :]}{1 / column_scale[j]}
        """

        n_rows, n_cols = X.shape

        column_means = np.asarray(column_means)
        if len(column_means) != n_cols:
            raise ValueError("Expected length %d but got shape %s" % (
                n_cols, column_means.shape))
        X = X - column_means.reshape((1, n_cols))
        column_weights = 1.0 / column_scales
        X *= column_weights.reshape((1, n_cols))
        row_means = np.zeros(n_rows, dtype=X.dtype)
        row_residual_sums = np.nansum(X, axis=1)
        for i in range(n_rows):
            row_mask = observed[i, :]
            sum_weights = column_weights[row_mask].sum()
            row_means[i] = row_residual_sums[i] / sum_weights
        return row_means