code: def truncated_normal_expval(mu, tau, a, b):
      
      phia = np.exp(normal_like(a, mu, tau))
      phib = np.exp(normal_like(b, mu, tau))
      sigma = 1. / np.sqrt(tau)
      Phia = utils.normcdf((a - mu) / sigma)
      if b == np.inf:
          Phib = 1.0
      else:
          Phib = utils.normcdf((b - mu) / sigma)
      return (mu + (phia - phib) / (Phib - Phia))[0]
docs: Expected value of the truncated normal distribution.
  
      .. math::
         E(X) =\mu + \frac{\sigma(\varphi_1-\varphi_2)}{T}
  
  
      where
  
      .. math::
         T & =\Phi\left(\frac{B-\mu}{\sigma}\right)-\Phi
         \left(\frac{A-\mu}{\sigma}\right)\text \\
         \varphi_1 &=
         \varphi\left(\frac{A-\mu}{\sigma}\right) \\
         \varphi_2 &=
         \varphi\left(\frac{B-\mu}{\sigma}\right) \\
  
      and :math:`\varphi = N(0,1)` and :math:`tau & 1/sigma**2`.
  
      :Parameters:
        - `mu` : Mean of the distribution.
        - `tau` : Precision of the distribution, which corresponds to 1/sigma**2 (tau > 0).
        - `a` : Left bound of the distribution.
        - `b` : Right bound of the distribution.