docs: def make_nn_funs(layer_sizes, L2_reg, noise_variance, nonlinearity=np.tanh):
      """These functions implement a standard multi-layer perceptron,
      vectorized over both training examples and weight samples."""
      shapes = list(zip(layer_sizes[:-1], layer_sizes[1:]))
      num_weights = sum((m+1)*n for m, n in shapes)
  
      def unpack_layers(weights):
          num_weight_sets = len(weights)
          for m, n in shapes:
              yield weights[:, :m*n]     .reshape((num_weight_sets, m, n)),\
                    weights[:, m*n:m*n+n].reshape((num_weight_sets, 1, n))
              weights = weights[:, (m+1)*n:]
  
      def predictions(weights, inputs):
          """weights is shape (num_weight_samples x num_weights)
             inputs  is shape (num_datapoints x D)"""
          inputs = np.expand_dims(inputs, 0)
          for W, b in unpack_layers(weights):
              outputs = np.einsum('mnd,mdo->mno', inputs, W) + b
              inputs = nonlinearity(outputs)
          return outputs
  
      def logprob(weights, inputs, targets):
          log_prior = -L2_reg * np.sum(weights**2, axis=1)
          preds = predictions(weights, inputs)
          log_lik = -np.sum((preds - targets)**2, axis=1)[:, 0] / noise_variance
          return log_prior + log_lik
  
      return num_weights, predictions, logprob