code: def get_lasso_coef_df_by_multi_sample_split(info_dict):
      """Gets the lasso coefficients' p-values and confidence intervals by multi sample-splitting.
  
      Repeats the single sample-splitting in
      `~greykite.algo.common.model_summary_utils.get_lasso_coef_by_single_sample_split`
      for multiple times, and aggregate the p-values and confidence intervals according to
      "High-Dimensional Inference: Confidence Intervals, p-Values and R-Software hdi"
      by Dezeure, Buhlmann, Meier and Meinshausen
  
      Parameters
      ----------
      info_dict : `dict`
          The information dictionary returned by
          `~greykite.algo.common.model_summary_utils.add_beta_var_lm`.
          Note that the "ml_model" in ``info_dict`` should be a lasso method
          that supports cloning.
  
      Returns
      -------
      coef_summary_df : `pandas.DataFrame`
          The coefficients summary dataframe.
      significance_code_legend : `str`
          The significance code legend.
      """
      # Keeping ``num_splits`` as 100 is a good choice.
      # Less number of splits make some values inaccurate (prob_nonzero, p_value, etc.).
      # Greater number of splits takes longer to run (each split runs a LassoCV).
      num_splits = 100
      beta_estimates = None
      standard_errs = None
      df_mses = np.array([])
      # beta_estimates includes the betas estimated by multi sample-splitting.
      # Each row is a coefficient, and each column is an estimation from a sample-splitting.
      for i in range(num_splits):
          (
              beta,
              standard_err,
              df_mse,
          ) = get_lasso_coef_by_single_sample_split(info_dict)
          if beta_estimates is None:
              beta_estimates = beta.reshape(-1, 1)
          else:
              beta_estimates = np.concatenate(
                  [beta_estimates, beta.reshape(-1, 1)], axis=1
              )
          if standard_errs is None:
              standard_errs = standard_err.reshape(-1, 1)
          else:
              standard_errs = np.concatenate(
                  [standard_errs, standard_err.reshape(-1, 1)], axis=1
              )
          df_mses = np.append(df_mses, df_mse)
      # Calculates nonzero probability for each feature coefficient
      beta_nonzero_prob = (
          np.count_nonzero(beta_estimates, axis=1)
          / beta_estimates.shape[1]
      )
  
      # Calculates individual p-values
      beta_estimates_p_values = np.where(
          beta_estimates != 0,
          2 * t.cdf(-np.abs(beta_estimates / standard_errs), df_mses),
          1,
      )  # p-value is 1.0 for zero coefficients
      # Aggregates indivudual p-values
      # Reference: P-Values for High-dimensional Regression
      # by Meinshausen, Meier, Buhlmann
      # Link: https://arxiv.org/pdf/0811.2177.pdf
      gamma_min = 0.05  # a lower bound for the p-value percentile, typically 0.05
      gammas = np.arange(gamma_min, 1, 0.01)
      qjs = None
      # For each percentile gamma, the qj is defined as the gamma_th percentile of p-values divided by gamma
      # Set to 1 if the value exceeds 1.
      for gamma in gammas:
          qj = (
              np.percentile(
                  beta_estimates_p_values, q=gamma * 100, axis=1
              )
              / gamma
          )
          qj = np.where(
              qj < 1, qj, 1
          )  # set to 1 if the value exceeds 1
          if qjs is None:
              qjs = qj.reshape(-1, 1)
          else:
              qjs = np.concatenate([qjs, qj.reshape(-1, 1)], axis=1)
      p_values = np.min(qjs, axis=1) * (1 - np.log(gamma_min))
      p_values = np.where(p_values < 1, p_values, 1)
      sig_codes = [
          "***"
          if p < 0.001
          else "**"
          if p < 0.01
          else "*"
          if p < 0.05
          else "."
          if p < 0.1
          else " "
          for p in p_values
      ]
  
      # Calculates p-value rankings, plus one because we want ranking to start with 1
      p_value_rankings = (
          beta_estimates_p_values.argsort(axis=1) + 1
      ) / beta_estimates_p_values.shape[1]
      # Calculates the confidence intervals
      # Reference: High-Dimensional Inference: Confidence Intervals, p-Values and R-Software hdi
      # by Dezeure, Buhlmann, Meier and Meinshausen
      alpha = 0.05
      # Divides by 2 because we need half of the area on each side
      quantile = (
          1 - alpha * p_value_rankings / (1 - np.log(gamma_min)) / 2
      )
      # Raveled dfs are applied along columns
      ordinary_confidence_interval_lb = (
          beta_estimates
          - t.ppf(quantile, df_mses.ravel()) * standard_errs
      )
      ordinary_confidence_interval_ub = (
          beta_estimates
          + t.ppf(quantile, df_mses.ravel()) * standard_errs
      )
  
      confidence_intervals = []
      for i in range(beta_estimates.shape[0]):
          if (beta_estimates[i] == 0).all():
              confidence_intervals.append([0.0, 0.0])
          else:
              check_index = (p_value_rankings[i] > gamma_min) & (
                  beta_estimates[i] != 0
              )
              if check_index.sum() == 0:
                  confidence_intervals.append([0.0, 0.0])
              else:
                  lb = np.min(
                      ordinary_confidence_interval_lb[i, check_index]
                  )
                  ub = np.max(
                      ordinary_confidence_interval_ub[i, check_index]
                  )
                  confidence_intervals.append([lb, ub])
  
      # Makes summary df
      coef_summary_df = pd.DataFrame(
          {
              "Pred_col": info_dict["pred_cols"],
              "Estimate": info_dict["beta"],
              "Pr(>)_split": p_values,
              "sig. code": sig_codes,
              "95%CI": confidence_intervals,
              "Prob_nonzero": beta_nonzero_prob,
          }
      )
      return (
          coef_summary_df,
          "0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      )
  
imports: import numpy
  import pandas
  import statsmodels
  import statsmodels.api
  from scipy.stats import f
  from scipy.stats import norm
  from scipy.stats import t
  from sklearn.base import clone
  from sklearn.linear_model import Ridge
  from sklearn.model_selection import train_test_split
  from greykite.algo.common.col_name_utils import INTERCEPT
  from greykite.algo.common.col_name_utils import simplify_pred_cols
label_desc: Use image node_id 13 for calling a global function with example usage: get_lasso_coef_df_by_multi_sample_split(info_dict) and returns: coef_summary_df,  str