code: def windowed_variance(
      x, low_indices=None, high_indices=None, axis=0, name=None
  ):
      """Windowed estimates of variance.
  
      Computes variances among data in the Tensor `x` along the given windows:
  
        result[i] = variance(x[low_indices[i]:high_indices[i]+1])
  
      accurately and efficiently.  To wit, if K is the size of
      `low_indices` and `high_indices`, and `N` is the size of `x` along
      the given `axis`, the computation takes O(K + N) work, O(log(N))
      depth (the length of the longest series of operations that are
      performed sequentially), and only uses O(1) TensorFlow kernel
      invocations.  The underlying algorithm is an adaptation of the
      streaming reduction for accurate variance computations given in [1].
  
      This function can be useful for assessing the behavior over time of
      trailing-window estimators from some iterative process, such as the
      last half of an MCMC chain.
  
      Suppose `x` has shape `Bx + [N] + E`, where the `Bx` component has
      rank `axis`, and `low_indices` and `high_indices` broadcast to shape
      `[M]`.  Then each element of `low_indices` and `high_indices`
      must be between 0 and N+1, and the shape of the output will be
      `Bx + [M] + E`.  Batch shape in the indices is not currently supported.
  
      The default windows are
      `[0, 1), [1, 2), [1, 3), [2, 4), [2, 5), ...`
      This corresponds to analyzing `x` as though it were streaming, for
      example successive states of an MCMC sampler, and we were interested
      in the variance of the last half of the data at each point.
  
      Args:
        x: A numeric `Tensor` holding `N` samples along the given `axis`,
          whose windowed variances are desired.
        low_indices: An integer `Tensor` defining the lower boundary
          (inclusive) of each window.  Default: elementwise half of
          `high_indices`.
        high_indices: An integer `Tensor` defining the upper boundary
          (exclusive) of each window.  Must be broadcast-compatible with
          `low_indices`.  Default: `tf.range(1, N+1)`, i.e., N windows
          that each end in the corresponding datum from `x` (inclusive)`.
        axis: Scalar `Tensor` designating the axis holding samples.  This
          is the axis of `x` along which we take windows, and therefore
          the axis that `low_indices` and `high_indices` index into.
          Other axes are treated in batch.  Default value: `0` (leftmost
          dimension).
        name: Python `str` name prefixed to Ops created by this function.
          Default value: `None` (i.e., `'windowed_variance'`).
  
      Returns:
        variances: A numeric `Tensor` holding the windowed variances of
          `x` along the `axis` dimension.
  
      #### References
      [1]: Philippe Pebay. Formulas for Robust, One-Pass Parallel Computation of
           Covariances and Arbitrary-Order Statistical Moments. _Technical Report
           SAND2008-6212_, 2008.
           https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
  
      """
      with tf.name_scope(name or "windowed_variance"):
          x = tf.convert_to_tensor(x)
          (
              low_indices,
              high_indices,
              low_counts,
              high_counts,
          ) = _prepare_window_args(x, low_indices, high_indices, axis)
  
          # We have a problem with indexing: the standard convention demands
          # the low index be inclusive, and the high index be exclusive.
          # However, tf.cumsum and cumulative_variance both include the ith
          # element in the ith result, so to implement the standard convention
          # we have to either invoke exclusive variants of the above, or
          # index off by one element.  Luckily, we can do the latter
          # without indexing off the beginning, because the value we fetch
          # when low_indices[i] == 0 or high_indices[i] == 0 is irrelevant,
          # because it gets multiplied by 0 later anyway.
          # Note that exclusive cumsum doesn't work either, because we
          # allow high_indices to take the value N+1 meaning "all the data".
          def index_for_cumulative(indices):
              return tf.maximum(indices - 1, 0)
  
          cum_sums = tf.cumsum(x, axis=axis)
          low_sums = tf.gather(
              cum_sums, index_for_cumulative(low_indices), axis=axis
          )
          high_sums = tf.gather(
              cum_sums, index_for_cumulative(high_indices), axis=axis
          )
          cum_variances = cumulative_variance(x, sample_axis=axis)
          low_variances = tf.gather(
              cum_variances,
              index_for_cumulative(low_indices),
              axis=axis,
          )
          high_variances = tf.gather(
              cum_variances,
              index_for_cumulative(high_indices),
              axis=axis,
          )
  
          # This formula is the binary accurate variance merge from [1],
          # adapted to subtract and batched across the indexed counts, sums,
          # and variances.
          # As a reminder, [1] shows the following for multisets A, B:
          #   var(A u B) = (|A|*var(A) + |B|*var(B) + correction) / |A u B|
          # where
          #   correction = (mean(A) - mean(B))**2 * |A| * |B| / |A u B|
          # For each high_indices[i] and low_indices[i], if we let
          # A be the multiset x[low_indices[i]:high_indices[i]] and B
          # be the multiset x[0:low_indices[i]], then
          #   var(A u B) = cum_variances[high_indices[i]] = high_variances[i]
          #   var(B) = cum_variances[low_indices[i]] = low_variances[i]
          # and the below solves for var(A).
          # This formula can also be read as implementing the above variance
          # computation by "unioning" A u B with a notional "negative B"
          # multiset.
          counts = high_counts - low_counts  # |A|
          discrepancies = (
              _safe_average(high_sums, high_counts)
              - _safe_average(low_sums, low_counts)
          ) ** 2  # (mean(A u B) - mean(B))**2
          adjustments = (
              high_counts * (-low_counts) / counts
          )  # |A u B| * -|B| / |A|
          residuals = (
              high_variances * high_counts
              - low_variances * low_counts
              + adjustments * discrepancies
          )
          return _safe_average(residuals, counts)
  
imports: import numpy
  import tensorflow.compat.v2
  from tensorflow_probability.python.internal import assert_util
  from tensorflow_probability.python.internal import distribution_util
  from tensorflow_probability.python.internal import dtype_util
  from tensorflow_probability.python.internal import prefer_static
  from tensorflow_probability.python.internal import tensorshape_util
  from tensorflow_probability.python.math.generic import reduce_logmeanexp
label_desc: Use image node_id 8 for calling a global function with example usage: windowed_variance(x, low_indices, high_indices, axis, name) and returns: _safe_average, tf