code: def average_weir_cockerham_fst(g, subpops, blen, max_allele=None):
      
  
      # calculate per-variant values
      a, b, c = weir_cockerham_fst(g, subpops, max_allele=max_allele)
  
      # calculate overall estimate
      a_sum = np.nansum(a)
      b_sum = np.nansum(b)
      c_sum = np.nansum(c)
      fst = a_sum / (a_sum + b_sum + c_sum)
  
      # compute the numerator and denominator within each block
      num_bsum = moving_statistic(a, statistic=np.nansum, size=blen)
      den_bsum = moving_statistic(a + b + c, statistic=np.nansum, size=blen)
  
      # calculate the statistic values in each block
      vb = num_bsum / den_bsum
  
      # estimate standard error
      _, se, vj = jackknife((num_bsum, den_bsum),
                            statistic=lambda n, d: np.sum(n) / np.sum(d))
  
      return fst, se, vb, vj
docs: Estimate average Fst and standard error using the block-jackknife.
  
      Parameters
      ----------
      g : array_like, int, shape (n_variants, n_samples, ploidy)
          Genotype array.
      subpops : sequence of sequences of ints
          Sample indices for each subpopulation.
      blen : int
          Block size (number of variants).
      max_allele : int, optional
          The highest allele index to consider.
  
      Returns
      -------
      fst : float
          Estimated value of the statistic using all data.
      se : float
          Estimated standard error.
      vb : ndarray, float, shape (n_blocks,)
          Value of the statistic in each block.
      vj : ndarray, float, shape (n_blocks,)
          Values of the statistic from block-jackknife resampling.