code: def GMM1_lpdf(
      samples, weights, mus, sigmas, low=None, high=None, q=None
  ):
      def print_verbose(s, x):
          return print(f"GMM1_lpdf:{s}", x)
  
      verbose = 0
      samples, weights, mus, sigmas = list(
          map(np.asarray, (samples, weights, mus, sigmas))
      )
      if samples.size == 0:
          return np.asarray([])
      if weights.ndim != 1:
          raise TypeError("need vector of weights", weights.shape)
      if mus.ndim != 1:
          raise TypeError("need vector of mus", mus.shape)
      if sigmas.ndim != 1:
          raise TypeError("need vector of sigmas", sigmas.shape)
      assert len(weights) == len(mus) == len(sigmas)
      _samples = samples
      samples = _samples.flatten()
  
      if verbose:
          print_verbose("samples", set(samples))
          print_verbose("weights", weights)
          print_verbose("mus", mus)
          print_verbose("sigmas", sigmas)
          print_verbose("low", low)
          print_verbose("high", high)
          print_verbose("q", q)
  
      if low is None and high is None:
          p_accept = 1
      else:
          p_accept = np.sum(
              weights
              * (
                  normal_cdf(high, mus, sigmas)
                  - normal_cdf(low, mus, sigmas)
              )
          )
  
      if q is None:
          dist = samples[:, None] - mus
          mahal = (dist / np.maximum(sigmas, EPS)) ** 2
          # mahal shape is (n_samples, n_components)
          Z = np.sqrt(2 * np.pi * sigmas**2)
          coef = weights / Z / p_accept
          rval = logsum_rows(-0.5 * mahal + np.log(coef))
      else:
          prob = np.zeros(samples.shape, dtype="float64")
          for w, mu, sigma in zip(weights, mus, sigmas):
              if high is None:
                  ubound = samples + q / 2
              else:
                  ubound = np.minimum(samples + q / 2, high)
              if low is None:
                  lbound = samples - q / 2
              else:
                  lbound = np.maximum(samples - q / 2, low)
              # -- two-stage addition is slightly more numerically accurate
              inc_amt = w * normal_cdf(ubound, mu, sigma)
              inc_amt -= w * normal_cdf(lbound, mu, sigma)
              prob += inc_amt
          rval = np.log(prob) - np.log(p_accept)
  
      if verbose:
          print_verbose("rval:", dict(list(zip(samples, rval))))
  
      rval.shape = _samples.shape
      return rval
  
imports: import logging
  import time
  import numpy
  from scipy.special import erf
  from .None import pyll
  from .pyll import scope
  from .pyll.stochastic import implicit_stochastic
  from .base import miscs_to_idxs_vals
  from .base import miscs_update_idxs_vals
  from .None import rand
label_desc: Use image node_id 6 for calling a global function with example usage: GMM1_lpdf(samples, weights, mus, sigmas, low, high, q) and returns: rval, print, np