column_center[j] =
        sum{i in observed[:, j]}{
            (1 / row_scale[i]) * (X[i, j]) - row_center[i])
        }
        ------------------------------------------------------------
        sum{i in observed[:, j]}{1 / row_scale[i]}

def estimate_column_means(
            self,
            X,
            observed,
            row_means,
            row_scales):
        """
        column_center[j] =
        sum{i in observed[:, j]}{
            (1 / row_scale[i]) * (X[i, j]) - row_center[i])
        }
        ------------------------------------------------------------
        sum{i in observed[:, j]}{1 / row_scale[i]}
        """
        n_rows, n_cols = X.shape
        row_means = np.asarray(row_means)

        if len(row_means) != n_rows:
            raise ValueError("Expected length %d but got shape %s" % (
                n_rows, row_means.shape))
        column_means = np.zeros(n_cols, dtype=X.dtype)

        X = X - row_means.reshape((n_rows, 1))
        row_weights = 1.0 / row_scales
        X *= row_weights.reshape((n_rows, 1))
        col_residual_sums = np.nansum(X, axis=0)
        for j in range(n_cols):
            col_mask = observed[:, j]
            sum_weights = row_weights[col_mask].sum()
            column_means[j] = col_residual_sums[j] / sum_weights
        return column_means