code: def estimate_column_means(
              self,
              X,
              observed,
              row_means,
              row_scales):
          
          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
docs: 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]}