code: def compute_pi_samples(
      n, h, states, sigma, alpha, theta, mean_y, seed=0, n_samples=200
  ):
      samples = np.full(
          (h, n_samples), fill_value=np.nan, dtype=np.float32
      )
      # states: level, meany, An, Bn, mu
      smoothed, _, A, B, _ = states[-1]
      np.random.seed(seed)
      for i in range(n, n + h):
          samples[i - n] = smoothed + (1 - 1 / theta) * (
              A * ((1 - alpha) ** i)
              + B * (1 - (1 - alpha) ** (i + 1)) / alpha
          )
          samples[i - n] += np.random.normal(
              scale=sigma, size=n_samples
          )
          smoothed = alpha * samples[i - n] + (1 - alpha) * smoothed
          mean_y = (i * mean_y + samples[i - n]) / (i + 1)
          B = (
              (i - 1) * B + 6 * (samples[i - n] - mean_y) / (i + 1)
          ) / (i + 2)
          A = mean_y - B * (i + 2) / 2
      return samples
  
imports: import math
  from typing import Tuple
  import numpy
  from numba import njit
  from scipy.stats import norm
  from statsmodels.tsa.seasonal import seasonal_decompose
  from statsmodels.tsa.stattools import acf
  from .ets import restrict_to_bounds, results
  from .utils import _seasonal_naive, _repeat_val_seas, CACHE, NOGIL
label_desc: Use image node_id 14 for calling a global function with example usage: compute_pi_samples(n, h, states, sigma, alpha, theta, mean_y, seed, n_samples) and returns: samples