code: def mvsdist(data):
      """
      'Frozen' distributions for mean, variance, and standard deviation of data.
  
      Parameters
      ----------
      data : array_like
          Input array. Converted to 1-D using ravel.
          Requires 2 or more data-points.
  
      Returns
      -------
      mdist : "frozen" distribution object
          Distribution object representing the mean of the data.
      vdist : "frozen" distribution object
          Distribution object representing the variance of the data.
      sdist : "frozen" distribution object
          Distribution object representing the standard deviation of the data.
  
      See Also
      --------
      bayes_mvs
  
      Notes
      -----
      The return values from ``bayes_mvs(data)`` is equivalent to
      ``tuple((x.mean(), x.interval(0.90)) for x in mvsdist(data))``.
  
      In other words, calling ``<dist>.mean()`` and ``<dist>.interval(0.90)``
      on the three distribution objects returned from this function will give
      the same results that are returned from `bayes_mvs`.
  
      References
      ----------
      T.E. Oliphant, "A Bayesian perspective on estimating mean, variance, and
      standard-deviation from data", https://scholarsarchive.byu.edu/facpub/278,
      2006.
  
      Examples
      --------
      >>> from scipy import stats
      >>> data = [6, 9, 12, 7, 8, 8, 13]
      >>> mean, var, std = stats.mvsdist(data)
  
      We now have frozen distribution objects "mean", "var" and "std" that we can
      examine:
  
      >>> mean.mean()
      9.0
      >>> mean.interval(0.95)
      (6.6120585482655692, 11.387941451734431)
      >>> mean.std()
      1.1952286093343936
  
      """
      x = ravel(data)
      n = len(x)
      if n < 2:
          raise ValueError("Need at least 2 data-points.")
      xbar = x.mean()
      C = x.var()
      if n > 1000:  # gaussian approximations for large n
          mdist = distributions.norm(loc=xbar, scale=math.sqrt(C / n))
          sdist = distributions.norm(
              loc=math.sqrt(C), scale=math.sqrt(C / (2.0 * n))
          )
          vdist = distributions.norm(
              loc=C, scale=math.sqrt(2.0 / n) * C
          )
      else:
          nm1 = n - 1
          fac = n * C / 2.0
          val = nm1 / 2.0
          mdist = distributions.t(
              nm1, loc=xbar, scale=math.sqrt(C / nm1)
          )
          sdist = distributions.gengamma(val, -2, scale=math.sqrt(fac))
          vdist = distributions.invgamma(val, scale=fac)
      return mdist, vdist, sdist
  
imports: from __future__ import annotations
  import math
  import warnings
  from collections import namedtuple
  import numpy
  from numpy import isscalar, r_, log, around, unique, asarray, zeros, arange, sort, amin, amax, atleast_1d, sqrt, array, compress, pi, exp, ravel, count_nonzero, sin, cos, arctan2, hypot
  from scipy import optimize
  from scipy import special
  from .None import statlib
  from .None import stats
  from .stats import find_repeats, _contains_nan, _normtest_finish
  from .contingency import chi2_contingency
  from .None import distributions
  from ._distn_infrastructure import rv_generic
  from ._hypotests import _get_wilcoxon_distr
label_desc: Use image node_id 2 for calling a global function with example usage: mvsdist(data) and returns: mdist,  vdist,  sdist