code: def compute_bootci(
      x,
      y=None,
      func=None,
      method="cper",
      paired=False,
      confidence=0.95,
      n_boot=2000,
      decimals=2,
      seed=None,
      return_dist=False,
  ):
      """Bootstrapped confidence intervals of univariate and bivariate functions.
  
      Parameters
      ----------
      x : 1D-array or list
          First sample. Required for both bivariate and univariate functions.
      y : 1D-array, list, or None
          Second sample. Required only for bivariate functions.
      func : str or custom function
          Function to compute the bootstrapped statistic. Accepted string values are:
  
          * ``'pearson'``: Pearson correlation (bivariate, paired x and y)
          * ``'spearman'``: Spearman correlation (bivariate, paired x and y)
          * ``'cohen'``: Cohen d effect size (bivariate, paired or unpaired x and y)
          * ``'hedges'``: Hedges g effect size (bivariate, paired or unpaired x and y)
          * ``'mean'``: Mean (univariate = only x)
          * ``'std'``: Standard deviation (univariate)
          * ``'var'``: Variance (univariate)
      method : str
          Method to compute the confidence intervals (see Notes):
  
          * ``'cper'``: Bias-corrected percentile method (default)
          * ``'norm'``: Normal approximation with bootstrapped bias and standard error
          * ``'per'``: Simple percentile
      paired : boolean
          Indicates whether x and y are paired or not. For example, for correlation functions or
          paired T-test, x and y are assumed to be paired. Pingouin will resample the pairs
          (x_i, y_i) when paired=True, and resample x and y separately when paired=False.
          If paired=True, x and y must have the same number of elements.
      confidence : float
          Confidence level (0.95 = 95%)
      n_boot : int
          Number of bootstrap iterations. The higher, the better, the slower.
      decimals : int
          Number of rounded decimals.
      seed : int or None
          Random seed for generating bootstrap samples.
      return_dist : boolean
          If True, return the confidence intervals and the bootstrapped distribution (e.g. for
          plotting purposes).
  
      Returns
      -------
      ci : array
          Bootstrapped confidence intervals.
  
      Notes
      -----
      Results have been tested against the
      `bootci <https://www.mathworks.com/help/stats/bootci.html>`_ Matlab function.
  
      Since version 1.7, SciPy also includes a built-in bootstrap function
      :py:func:`scipy.stats.bootstrap`. The SciPy implementation has two advantages over Pingouin: it
      is faster when using ``vectorized=True``, and it supports the bias-corrected and accelerated
      (BCa) confidence intervals for univariate functions. However, unlike Pingouin, it does not
      return the bootstrap distribution.
  
      The percentile bootstrap method (``per``) is defined as the
      :math:`100 \\times \\frac{\\alpha}{2}` and :math:`100 \\times \\frac{1 - \\alpha}{2}`
      percentiles of the distribution of :math:`\\theta` estimates obtained from resampling, where
      :math:`\\alpha` is the level of significance (1 - confidence, default = 0.05 for 95% CIs).
  
      The bias-corrected percentile method (``cper``) corrects for bias of the bootstrap
      distribution. This method is different from the BCa method — the default in Matlab and SciPy —
      which corrects for both bias and skewness of the bootstrap distribution using jackknife
      resampling.
  
      The normal approximation method (``norm``) calculates the confidence intervals with the
      standard normal distribution using bootstrapped bias and standard error.
  
      References
      ----------
      * DiCiccio, T. J., & Efron, B. (1996). Bootstrap confidence intervals. Statistical science,
        189-212.
  
      * Davison, A. C., & Hinkley, D. V. (1997). Bootstrap methods and their application (Vol. 1).
        Cambridge university press.
  
      * Jung, Lee, Gupta, & Cho (2019). Comparison of bootstrap confidence interval methods for
        GSCA using a Monte Carlo simulation. Frontiers in psychology, 10, 2215.
  
      Examples
      --------
      1. Bootstrapped 95% confidence interval of a Pearson correlation
  
      >>> import pingouin as pg
      >>> import numpy as np
      >>> rng = np.random.default_rng(42)
      >>> x = rng.normal(loc=4, scale=2, size=100)
      >>> y = rng.normal(loc=3, scale=1, size=100)
      >>> stat = np.corrcoef(x, y)[0][1]
      >>> ci = pg.compute_bootci(x, y, func='pearson', paired=True, seed=42, decimals=4)
      >>> print(round(stat, 4), ci)
      0.0945 [-0.098   0.2738]
  
      Let's compare to SciPy's built-in bootstrap function
  
      >>> from scipy.stats import bootstrap
      >>> bt_scipy = bootstrap(
      ...       data=(x, y), statistic=lambda x, y: np.corrcoef(x, y)[0][1],
      ...       method="basic", vectorized=False, n_resamples=2000, paired=True, random_state=42)
      >>> np.round(bt_scipy.confidence_interval, 4)
      array([-0.0952,  0.2883])
  
      2. Bootstrapped 95% confidence interval of a Cohen d
  
      >>> stat = pg.compute_effsize(x, y, eftype='cohen')
      >>> ci = pg.compute_bootci(x, y, func='cohen', seed=42, decimals=3)
      >>> print(round(stat, 4), ci)
      0.7009 [0.403 1.009]
  
      3. Bootstrapped confidence interval of a standard deviation (univariate)
  
      >>> import numpy as np
      >>> stat = np.std(x, ddof=1)
      >>> ci = pg.compute_bootci(x, func='std', seed=123)
      >>> print(round(stat, 4), ci)
      1.5534 [1.38 1.8 ]
  
      Compare to SciPy's built-in bootstrap function, which returns the bias-corrected and
      accelerated CIs (see Notes).
  
      >>> def std(x, axis):
      ...     return np.std(x, ddof=1, axis=axis)
      >>> bt_scipy = bootstrap(data=(x, ), statistic=std, n_resamples=2000, random_state=123)
      >>> np.round(bt_scipy.confidence_interval, 2)
      array([1.39, 1.81])
  
      Changing the confidence intervals type in Pingouin
  
      >>> pg.compute_bootci(x, func='std', seed=123, method="norm")
      array([1.37, 1.76])
  
      >>> pg.compute_bootci(x, func='std', seed=123, method="percentile")
      array([1.35, 1.75])
  
      4. Bootstrapped confidence interval using a custom univariate function
  
      >>> from scipy.stats import skew
      >>> round(skew(x), 4), pg.compute_bootci(x, func=skew, n_boot=10000, seed=123)
      (-0.137, array([-0.55,  0.32]))
  
      5. Bootstrapped confidence interval using a custom bivariate function. Here, x and y are not
      paired and can therefore have different sizes.
  
      >>> def mean_diff(x, y):
      ...     return np.mean(x) - np.mean(y)
      >>> y2 = rng.normal(loc=3, scale=1, size=200)  # y2 has 200 samples, x has 100
      >>> ci = pg.compute_bootci(x, y2, func=mean_diff, n_boot=10000, seed=123)
      >>> print(round(mean_diff(x, y2), 2), ci)
      0.88 [0.54 1.21]
  
      We can also get the bootstrapped distribution
  
      >>> ci, bt = pg.compute_bootci(x, y2, func=mean_diff, n_boot=10000, return_dist=True, seed=9)
      >>> print(f"The bootstrap distribution has {bt.size} samples. The mean and standard "
      ...       f"{bt.mean():.4f} ± {bt.std():.4f}")
      The bootstrap distribution has 10000 samples. The mean and standard 0.8807 ± 0.1704
      """
      from inspect import isfunction, isroutine
      from scipy.stats import norm
  
      # Check other arguments
      assert isinstance(confidence, float)
      assert 0 < confidence < 1, "confidence must be between 0 and 1."
      assert method in [
          "norm",
          "normal",
          "percentile",
          "per",
          "cpercentile",
          "cper",
      ]
      assert (
          isfunction(func) or isinstance(func, str) or isroutine(func)
      ), (
          "func must be a function (e.g. np.mean, custom function) or a string (e.g. 'pearson'). "
          "See documentation for more details."
      )
      vectorizable = False
  
      # Check x
      x = np.asarray(x)
      nx = x.size
      assert x.ndim == 1, "x must be one-dimensional."
      assert nx > 1, "x must have more than one element."
  
      # Check y
      if y is not None:
          y = np.asarray(y)
          ny = y.size
          assert y.ndim == 1, "y must be one-dimensional."
          assert ny > 1, "y must have more than one element."
          if paired:
              assert (
                  nx == ny
              ), "x and y must have the same number of elements when paired=True."
  
      # Check string functions
      if isinstance(func, str):
          func_str = "%s" % func
          if func == "pearson":
              assert (
                  paired
              ), "Paired should be True if using correlation functions."
  
              def func(x, y):
                  return pearsonr(x, y)[0]  # Faster than np.corrcoef
  
          elif func == "spearman":
              from scipy.stats import spearmanr
  
              assert (
                  paired
              ), "Paired should be True if using correlation functions."
  
              def func(x, y):
                  return spearmanr(x, y)[0]
  
          elif func in ["cohen", "hedges"]:
              from pingouin.effsize import compute_effsize
  
              def func(x, y):
                  return compute_effsize(
                      x, y, paired=paired, eftype=func_str
                  )
  
          elif func == "mean":
              vectorizable = True
  
              def func(x):
                  return np.mean(x, axis=0)
  
          elif func == "std":
              vectorizable = True
  
              def func(x):
                  return np.std(x, ddof=1, axis=0)
  
          elif func == "var":
              vectorizable = True
  
              def func(x):
                  return np.var(x, ddof=1, axis=0)
  
          else:
              raise ValueError("Function string not recognized.")
  
      # Bootstrap
      bootstat = np.empty(n_boot)
      rng = np.random.default_rng(seed)  # Random seed
      boot_x = rng.choice(
          np.arange(nx), size=(nx, n_boot), replace=True
      )
  
      if y is not None:
          reference = func(x, y)
          if paired:
              for i in range(n_boot):
                  # Note that here we use a bootstrapping procedure with replacement
                  # of all the pairs (Xi, Yi). This is NOT suited for
                  # hypothesis testing such as p-value estimation). Instead, for the
                  # latter, one must only shuffle the Y values while keeping the X
                  # values constant, i.e.:
                  # >>> boot_x = rng.random_sample((n_boot, n)).argsort(axis=1)
                  # >>> for i in range(n_boot):
                  # >>>   bootstat[i] = func(x, y[boot_x[i, :]])
                  bootstat[i] = func(x[boot_x[:, i]], y[boot_x[:, i]])
          else:
              boot_y = rng.choice(
                  np.arange(ny), size=(ny, n_boot), replace=True
              )
              for i in range(n_boot):
                  bootstat[i] = func(x[boot_x[:, i]], y[boot_y[:, i]])
      else:
          reference = func(x)
          if vectorizable:
              bootstat = func(x[boot_x])
          else:
              for i in range(n_boot):
                  bootstat[i] = func(x[boot_x[:, i]])
  
      # CONFIDENCE INTERVALS
      # See Matlab bootci function
      alpha = (1 - confidence) / 2
      if method in ["norm", "normal"]:
          # Normal approximation
          za = norm.ppf(alpha)  # = 1.96
          se = np.std(bootstat, ddof=1)
          bias = np.mean(bootstat - reference)
          ci = np.array(
              [reference - bias + se * za, reference - bias - se * za]
          )
      elif method in ["percentile", "per"]:
          # Simple percentile
          interval = 100 * np.array([alpha, 1 - alpha])
          ci = np.percentile(bootstat, interval)
          pass
      else:
          # Bias-corrected percentile bootstrap
          from pingouin.regression import _bias_corrected_ci
  
          ci = _bias_corrected_ci(
              bootstat, reference, alpha=(1 - confidence)
          )
  
      ci = np.round(ci, decimals)
  
      if return_dist:
          return ci, bootstat
      else:
          return ci
  
file_path: pingouin/pingouin/effsize.py
imports: import warnings
  import numpy
  from scipy.stats import pearsonr
  from pingouin.utils import _check_eftype, remove_na
label_desc: Use image node_id 2 for calling a global function with example usage: compute_bootci(x, y, func, method, paired, confidence, n_boot, decimals, seed, return_dist) and returns: ci,  bootstat, ci, pearsonr, spearmanr, compute_effsize, np, np, np