code: def bvnmom(
      r: Tensor,
      xl: Tensor,
      yl: Tensor,
      xu: Tensor,
      yu: Tensor,
      p: Optional[Tensor] = None,
  ) -> Tuple[Tensor, Tensor]:
      r"""Computes the expected values of truncated, bivariate normal random variables.
  
      Let `x` and `y` be a pair of standard bivariate normal random variables having
      correlation `r`. This function computes `E([x,y] \| [xl,yl] < [x,y] < [xu,yu])`.
  
      Following [Muthen1990moments]_ equations (4) and (5), we have
  
          `E(x \| [xl, yl] < [x, y] < [xu, yu]) \
          = Z^{-1} \phi(xl) P(yl < y < yu \| x=xl) - \phi(xu) P(yl < y < yu \| x=xu),`
  
      where `Z = P([xl, yl] < [x, y] < [xu, yu])` and `\phi` is the standard normal PDF.
  
      Args:
          r: Tensor of correlation coefficients.
          xl: Tensor of lower bounds for `x`, same shape as `r`.
          xu: Tensor of upper bounds for `x`, same shape as `r`.
          yl: Tensor of lower bounds for `y`, same shape as `r`.
          yu: Tensor of upper bounds for `y`, same shape as `r`.
          p: Tensor of probabilities `P(xl < x < xu, yl < y < yu)`, same shape as `r`.
  
      Returns:
          `E(x \| [xl, yl] < [x, y] < [xu, yu])` and
          `E(y \| [xl, yl] < [x, y] < [xu, yu])`.
      """
      if not (r.shape == xl.shape == xu.shape == yl.shape == yu.shape):
          raise UnsupportedError(
              "Arguments to `bvn` must have the same shape."
          )
  
      if p is None:
          p = bvn(r=r, xl=xl, xu=xu, yl=yl, yu=yu)
  
      corr = r[..., None, None]
      istd = (1 - corr.square()).rsqrt()
      lower = torch.stack([xl, yl], -1)
      upper = torch.stack([xu, yu], -1)
      bounds = torch.stack([lower, upper], -1)
      deltas = safe_mul(corr, bounds)
  
      # Compute densities and conditional probabilities
      density_at_bounds = phi(bounds)
      prob_given_bounds = Phi(
          safe_mul(istd, safe_sub(upper.flip(-1).unsqueeze(-1), deltas))
      ) - Phi(
          safe_mul(istd, safe_sub(lower.flip(-1).unsqueeze(-1), deltas))
      )
  
      # Evaluate Muthen's formula
      p_diffs = (
          -(density_at_bounds * prob_given_bounds).diff().squeeze(-1)
      )
      moments = (1 / p).unsqueeze(-1) * (
          p_diffs + r.unsqueeze(-1) * p_diffs.flip(-1)
      )
      return moments.unbind(-1)
  
file_path: botorch/botorch/utils/probability/bvn.py
imports: from __future__ import annotations
  from math import pi
  from typing import Optional, Tuple
  import torch
  from botorch.exceptions import UnsupportedError
  from botorch.utils.probability.utils import case_dispatcher, get_constants_like, leggauss, ndtr, phi, STANDARDIZED_RANGE
  from botorch.utils.safe_math import div, exp, mul, sub
  from torch import Tensor
label_desc: Use image node_id 5 for calling a global function with example usage: bvnmom(r, xl, yl, xu, yu, p) and returns: moments