code: def point_predict(f, x, size=1, nugget=None):
      
      orig_shape = x.shape
      x = regularize_array(x)
  
      mu = f.M_internal(x, regularize=False)
      V = f.C_internal(x, regularize=False)
      if nugget is not None:
          V += nugget
      out = random.normal(size=(size, x.shape[0])) * sqrt(V) + mu
      return out.reshape((size,) + orig_shape[:-1]).squeeze()
docs: point_predict(f, x[, size, nugget])
  
      Makes 'size' simulations for f(x) + N(0,nugget).
      Simulated values of f(x_i) are uncorrelated for different i.
      Useful for geostatistical predictions.