### Hyper-Parameter Optimization Output Example

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/MVN_LowRank.ipynb

Example output from the hyper-parameter optimization process, showing study creation, progress, and trial results with optimized parameters and their corresponding values.

```text
[I 2023-06-21 16:08:24,386] A new study created in memory with name: XGBoostLSS Hyper-Parameter Optimization\nC:\Users\maerzale\.virtualenvs\XGBoostLSS-vIPRRz-M\lib\site-packages\optuna\progress_bar.py:56: ExperimentalWarning: Progress bar is experimental (supported from v1.2.0). The interface can change in the future.\n\n  0%|          | 0/20 [00:00<?, ?it/s]\n\n[I 2023-06-21 16:10:58,756] Trial 0 finished with value: 4715.1293944 and parameters: {'eta': 0.011583802021120177, 'max_depth': 6, 'gamma': 4.298718115769394e-07, 'subsample': 0.4389977358989794, 'colsample_bytree': 0.2937778261408931, 'min_child_weight': 2.0751720031622087e-06, 'booster': 'gbtree'}. Best is trial 0 with value: 4715.1293944.\n[I 2023-06-21 16:13:46,045] Trial 1 finished with value: 5506.9937502 and parameters: {'eta': 0.0006950905309493148, 'max_depth': 9, 'gamma': 0.0005319868904179334, 'subsample': 0.22711708632866234, 'colsample_bytree': 0.6606216557241239, 'min_child_weight': 1.7073883542887016e-07, 'booster': 'gbtree'}. Best is trial 0 with value: 4715.1293944.\n[I 2023-06-21 16:16:39,844] Trial 2 finished with value: 5601.6347656 and parameters: {'eta': 0.00010131220313381806, 'max_depth': 9, 'gamma': 6.680066899855348e-06, 'subsample': 0.8351289737147574, 'colsample_bytree': 0.6594140272475766, 'min_child_weight': 1.990513192590756e-06, 'booster': 'gbtree'}. Best is trial 0 with value: 4715.1293944.\n[I 2023-06-21 16:19:11,556] Trial 3 finished with value: 5299.090625 and parameters: {'eta': 0.0022487377373350856, 'max_depth': 4, 'gamma': 0.002434513612715828, 'subsample': 0.4716580160625341, 'colsample_bytree': 0.30539812876010936, 'min_child_weight': 23.460455962477635, 'booster': 'gbtree'}. Best is trial 0 with value: 4715.1293944.
```

--------------------------------

### Install XGBoostLSS Development Version

Source: https://github.com/statmixedml/xgboostlss/blob/master/README.md

Use this command to install the latest development version of XGBoostLSS directly from GitHub.

```bash
pip install git+https://github.com/StatMixedML/XGBoostLSS.git
```

--------------------------------

### Install XGBoostLSS PyPI Version

Source: https://github.com/statmixedml/xgboostlss/blob/master/README.md

Use this command to install the stable version of XGBoostLSS from the Python Package Index (PyPI).

```bash
pip install xgboostlss
```

--------------------------------

### Install SHAP dependency

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/index.md

Install the SHAP dependency required for model explanation features.

```python
pip install git+https://github.com/dsgibbons/shap.git
```

--------------------------------

### Hyperparameter Optimization with hyper_opt()

Source: https://context7.com/statmixedml/xgboostlss/llms.txt

Automates hyperparameter tuning using Optuna. Requires the xgboostlss[all_extras] installation.

```python
from xgboostlss.model import XGBoostLSS
from xgboostlss.distributions.Gaussian import Gaussian
from xgboostlss.datasets.data_loader import load_simulated_gaussian_data
import xgboost as xgb
import numpy as np

# Load data
train, _ = load_simulated_gaussian_data()
dtrain = xgb.DMatrix(train.filter(regex="x"), label=train["y"].values)

# Initialize model
xgblss = XGBoostLSS(Gaussian())

# Define hyperparameter search space
param_dict = {
    "eta":              ["float", {"low": 1e-5,   "high": 1,     "log": True}],
    "max_depth":        ["int",   {"low": 1,      "high": 10,    "log": False}],
    "gamma":            ["float", {"low": 1e-8,   "high": 40,    "log": True}],
    "subsample":        ["float", {"low": 0.2,    "high": 1.0,   "log": False}],
    "colsample_bytree": ["float", {"low": 0.2,    "high": 1.0,   "log": False}],
    "min_child_weight": ["float", {"low": 1e-8,   "high": 500,   "log": True}],
    "booster":          ["categorical", ["gbtree"]]
}

# Run hyperparameter optimization
np.random.seed(123)
opt_params = xgblss.hyper_opt(
    hp_dict=param_dict,
    dtrain=dtrain,
    num_boost_round=100,      # Max boosting iterations per trial
    nfold=5,                  # CV folds
    early_stopping_rounds=20, # Early stopping rounds
    max_minutes=10,           # Time budget
    n_trials=30,              # Number of trials
    silence=False,            # Show progress
    seed=123,                 # CV seed
    hp_seed=123               # Bayesian search seed
)

# opt_params contains best parameters including optimal number of rounds
print(f"Best parameters: {opt_params}")
# Output: {'eta': 0.09, 'max_depth': 2, ..., 'opt_rounds': 52}

# Train final model with optimal parameters
n_rounds = opt_params.pop("opt_rounds")
xgblss.train(opt_params, dtrain, num_boost_round=n_rounds)
```

--------------------------------

### Optimization Output Logs

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/SplineFlow_Regression.ipynb

Example console output generated during the hyperparameter optimization process.

```text
[I 2023-07-20 15:21:25,613] A new study created in memory with name: XGBoostLSS Hyper-Parameter Optimization
```

```text
  0%|          | 0/100 [00:00<?, ?it/s]
```

```text
[I 2023-07-20 15:22:30,946] Trial 0 finished with value: 6255.344043 and parameters: {'eta': 0.03692566115300059, 'max_depth': 4, 'gamma': 7.887574463976319e-08, 'subsample': 0.9144017805889704, 'colsample_bytree': 0.29974692606946474, 'min_child_weight': 2.1274056292372507e-07, 'booster': 'gbtree'}. Best is trial 0 with value: 6255.344043.
[I 2023-07-20 15:22:44,079] Trial 1 finished with value: 16614.029491999998 and parameters: {'eta': 0.0017154637091307261, 'max_depth': 7, 'gamma': 8.23346450464663e-07, 'subsample': 0.5893614292128422, 'colsample_bytree': 0.22046879379684717, 'min_child_weight': 1.0217730512063771e-06, 'booster': 'gbtree'}. Best is trial 0 with value: 6255.344043.
```

--------------------------------

### Import Libraries for XGBoostLSS

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/MVN_LowRank.ipynb

Imports necessary libraries for XGBoostLSS, data loading, manipulation, and plotting. Ensure these are installed before running the code.

```python
from xgboostlss.model import *
from xgboostlss.distributions.MVN_LoRa import *
from xgboostlss.datasets.data_loader import load_simulated_multivariate_gaussian_data

from sklearn.model_selection import train_test_split
import pandas as pd
import multiprocessing
import plotnine
from plotnine import *
plotnine.options.figure_size = (18, 9)
n_cpu = multiprocessing.cpu_count()
```

--------------------------------

### Import necessary libraries for XGBoostLSS and data handling

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/GaussianMixture_Regression_CaliforniaHousing.ipynb

Imports core components from xgboostlss, scikit-learn for data, and libraries for multiprocessing and plotting. Ensure these libraries are installed before running.

```python
from xgboostlss.model import *
from xgboostlss.distributions.Gaussian import *
from xgboostlss.distributions.Mixture import *
from xgboostlss.distributions.mixture_distribution_utils import MixtureDistributionClass

from sklearn import datasets
from sklearn.model_selection import train_test_split
import multiprocessing
import matplotlib.pyplot as plt
import seaborn as sns
figure_size = (10,5)
import plotnine
from plotnine import *
plotnine.options.figure_size = figure_size
```

--------------------------------

### Define Hyper-Parameter Search Space

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/ZAGamma_Regression.ipynb

Defines a dictionary specifying the parameters to tune, their sampling types (float, int, categorical, none), and their ranges or choices. This setup is crucial for guiding the hyper-parameter optimization process.

```python
param_dict = {
    "eta":              ["float", {"low": 1e-5,   "high": 1,     "log": True}],
    "max_depth":        ["int",   {"low": 1,      "high": 10,    "log": False}],
    "gamma":            ["float", {"low": 1e-8,   "high": 40,    "log": True}],
    "subsample":        ["float", {"low": 0.2,    "high": 1.0,   "log": False}],
    "colsample_bytree": ["float", {"low": 0.2,    "high": 1.0,   "log": False}],
    "booster":          ["categorical", ["gbtree"]],
    # "tree_method":    ["categorical", ["auto", "approx", "hist", "gpu_hist"]],
    # "gpu_id":         ["none", [0]]
}
```

--------------------------------

### Optimization Process Logs

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Dirichlet_Regression.ipynb

Console output showing the initialization and progress of the Optuna study.

```text
[I 2023-06-21 16:12:29,791] A new study created in memory with name: XGBoostLSS Hyper-Parameter Optimization
C:\Users\maerzale\.virtualenvs\XGBoostLSS-vIPRRz-M\lib\site-packages\optuna\progress_bar.py:56: ExperimentalWarning: Progress bar is experimental (supported from v1.2.0). The interface can change in the future.
```

```text
  0%|          | 0/20 [00:00<?, ?it/s]
```

```text
[I 2023-06-21 16:12:32,258] Trial 0 finished with value: -19.765962000000002 and parameters: {'eta': 3.5623752946431085e-07, 'max_depth': 1, 'gamma': 4.0777764217919787e-08, 'min_child_weight': 3.2270050801713505e-06, 'booster': 'gbtree'}. Best is trial 0 with value: -19.765962000000002.
[I 2023-06-21 16:12:33,572] Trial 1 finished with value: -39.1656265 and parameters: {'eta': 0.08832668402281392, 'max_depth': 2, 'gamma': 2.103048813930706e-05, 'min_child_weight': 0.0002905776860050202, 'booster': 'gbtree'}. Best is trial 1 with value: -39.1656265.
[I 2023-06-21 16:12:35,611] Trial 2 finished with value: -19.772727 and parameters: {'eta': 2.2096866886454325e-06, 'max_depth': 3, 'gamma': 1.8295275180845298, 'min_child_weight': 1.7565392988395257e-08, 'booster': 'gbtree'}. Best is trial 1 with value: -39.1656265.
[I 2023-06-21 16:12:37,404] Trial 3 finished with value: -19.764789999999998 and parameters: {'eta': 5.080034826525112e-08, 'max_depth': 4, 'gamma': 8.693694709232917, 'min_child_weight': 0.04959110160428548, 'booster': 'gbtree'}. Best is trial 1 with value: -39.1656265.
[I 2023-06-21 16:12:39,156] Trial 4 finished with value: -19.764696 and parameters: {'eta': 1.1568572916751665e-08, 'max_depth': 3, 'gamma': 0.0004128157890043689, 'min_child_weight': 5.372767848252571, 'booster': 'gbtree'}. Best is trial 1 with value: -39.1656265.
```

--------------------------------

### Initialize XGBoostLSS Environment

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/ZAGamma_Regression.ipynb

Imports necessary libraries and configures the environment for XGBoostLSS modeling.

```python
from xgboostlss.model import *
from xgboostlss.distributions.ZAGamma import *

from sklearn.model_selection import train_test_split 
import pandas as pd
import multiprocessing
import plotnine
from plotnine import *
plotnine.options.figure_size = (18, 9)
n_cpu = multiprocessing.cpu_count()
```

--------------------------------

### Importing XGBoostLSS and dependencies

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Gaussian_Regression.ipynb

Initializes the environment by importing model classes, distribution definitions, and data loading utilities.

```python
from xgboostlss.model import *
from xgboostlss.distributions.Gaussian import *
from xgboostlss.datasets.data_loader import load_simulated_gaussian_data
from scipy.stats import norm

import multiprocessing
import plotnine
from plotnine import *
plotnine.options.figure_size = (12, 8)
```

--------------------------------

### Initialize XGBoostLSS with Expectile Distribution

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Expectile_Regression.ipynb

Initializes the XGBoostLSS model specifying the Expectile distribution. Key parameters include stabilization method, the list of expectiles to estimate, and whether to penalize crossing of expectiles.

```python
xgblss = XGBoostLSS(
    Expectile(stabilization="None",              # Options are "None", "MAD", "L2".
              expectiles = [0.05, 0.95],         # List of expectiles to be estimated, in increasing order.
              penalize_crossing = True           # Whether to include a penalty term to discourage crossing of expectiles.
              )    
)
```

--------------------------------

### Configuring the Gaussian distribution model

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Gaussian_Regression.ipynb

Initializes the XGBoostLSS model with specific stabilization, response function, and loss function settings.

```python
xgblss = XGBoostLSS(
    Gaussian(stabilization="None",  
             response_fn="exp",      
             loss_fn="nll"          
            )
)
```

--------------------------------

### Predict Distributional Parameters with XGBoostlss

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Dirichlet_Regression.ipynb

Use the `predict` method with `pred_type='parameters'` to get the raw distributional parameters. These parameters are not normalized to the simplex by default.

```python
pred_params = xgblss.predict(dtrain,
                             pred_type="parameters")
```

--------------------------------

### Get Quantile Predictions from XGBoostlss

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/SplineFlow_Regression.ipynb

Use the predict method with pred_type='quantiles' to obtain predicted quantiles for specified samples. Ensure n_samples and quantiles are set appropriately.

```python
pred_quantiles = xgblss.predict(dtest,
                                pred_type="quantiles",
                                n_samples=n_samples,
                                quantiles=quant_sel)
```

--------------------------------

### Initialize XGBoostLSS Distribution

Source: https://context7.com/statmixedml/xgboostlss/llms.txt

Configures the distribution object with stabilization, response, and loss function parameters.

```python
dist = Gaussian(
    stabilization="None",  # "None", "MAD", "L2" - gradient stabilization
    response_fn="exp",     # "exp", "softplus" - ensures positive scale
    loss_fn="nll",         # "nll", "crps" - loss function
    initialize=False       # Use MLE starting values
)
```

--------------------------------

### Initialize XGBoostLSS with MVN_LoRa Distribution

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/MVN_LowRank.ipynb

Initializes an XGBoostLSS model using the Multivariate Normal distribution with Low-Rank Approximation (MVN_LoRa). Key parameters include the number of targets (D), the rank of the approximation, stabilization method, response function, and loss function.

```python
# Specifies a multivariate Normal distribution, using the Low-Rank approximation. See ?MVN_LoRa for details.
xgblss = XGBoostLSS(
    MVN_LoRa(D=n_targets,                      # Specifies the number of targets
             rank=2,                           # Specifies the rank 
             stabilization="None",             # Options are "None", "MAD", "L2".
             response_fn="exp",                # Function to transform the lower-triangular factor of the covariance, e.g., "exp" or "softplus".
             loss_fn="nll"                     # Loss function, i.e., nll.
            )               
)
```

--------------------------------

### Plot SHAP Feature Importance and Partial Dependence

Source: https://context7.com/statmixedml/xgboostlss/llms.txt

Use this to visualize feature importance and partial dependence for specific distributional parameters. Ensure the 'shap' library is installed.

```python
from xgboostlss.model import XGBoostLSS
from xgboostlss.distributions.Gaussian import Gaussian
from xgboostlss.datasets.data_loader import load_simulated_gaussian_data
import xgboost as xgb

# Load and train model
train, test = load_simulated_gaussian_data()
dtrain = xgb.DMatrix(train.filter(regex="x"), label=train["y"].values)
X_test = test.filter(regex="x")

xgblss = XGBoostLSS(Gaussian())
xgblss.train({"eta": 0.1, "max_depth": 3}, dtrain, num_boost_round=50)

# Plot partial dependence of scale parameter on x_true feature
xgblss.plot(
    X=X_test,
    feature="x_true",           # Feature to analyze
    parameter="scale",          # Distributional parameter (loc, scale)
    plot_type="Partial_Dependence"
)

# Plot feature importance for scale parameter
xgblss.plot(
    X=X_test,
    parameter="scale",
    max_display=15,             # Show top 15 features
    plot_type="Feature_Importance"
)

# For location parameter
xgblss.plot(X=X_test, parameter="loc", plot_type="Feature_Importance")
```

--------------------------------

### Import Libraries and Load Data

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Expectile_Regression.ipynb

Imports necessary libraries and loads simulated Gaussian data. It then preprocesses the data into training and testing sets, creating DMatrix objects for XGBoost.

```python
from xgboostlss.model import *
from xgboostlss.distributions.Expectile import *
from xgboostlss.datasets.data_loader import load_simulated_gaussian_data

import multiprocessing
import plotnine
from plotnine import *
plotnine.options.figure_size = (20, 10)
```

```python
# The data is a simulated Gaussian as follows, where x is the only true feature and all others are noise variables
    # loc = 10
    # scale = 1 + 4*((0.3 < x) & (x < 0.5)) + 2*(x > 0.7)

train, test = load_simulated_gaussian_data()
n_cpu = multiprocessing.cpu_count()

X_train, y_train = train.filter(regex="x"), train["y"].values
X_test, y_test = test.filter(regex="x"), test["y"].values

dtrain = xgb.DMatrix(X_train, label=y_train, nthread=n_cpu)
dtest = xgb.DMatrix(X_test, nthread=n_cpu)
```

--------------------------------

### Configure and Execute Hyper-Parameter Optimization

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Dirichlet_Regression.ipynb

Defines the parameter search space and triggers the optimization process using xgblss.hyper_opt.

```python
param_dict = {
    "eta":              ["float", {"low": 1e-08,  "high": 1,     "log": True}],
    "max_depth":        ["int",   {"low": 1,      "high": 5,     "log": False}],
    "gamma":            ["float", {"low": 1e-8,   "high": 40,    "log": True}],
    "min_child_weight": ["float", {"low": 1e-8,   "high": 500,   "log": True}]
}

np.random.seed(123)
opt_param = xgblss.hyper_opt(param_dict,
                             dtrain,
                             num_boost_round=100,        # Number of boosting iterations.
                             nfold=2,                    # Number of cv-folds.
                             early_stopping_rounds=20,   # Number of early-stopping rounds
                             max_minutes=120,            # Time budget in minutes, i.e., stop study after the given number of minutes.
                             n_trials=20,                # The number of trials. If this argument is set to None, there is no limitation on the number of trials.
                             silence=False,              # Controls the verbosity of the trail, i.e., user can silence the outputs of the trail.
                             seed=123,                   # Seed used to generate cv-folds.
                             hp_seed=None                # Seed for random number generator used in the Bayesian hyperparameter search.
                            )
```

--------------------------------

### Initialize XGBoostLSS with Spline-Flow

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/SplineFlow_Regression.ipynb

Configures an XGBoostLSS model using the Spline-Flow distribution with specific target support and spline parameters.

```python
bound = np.max([np.abs(y_train.min()), y_train.max()])

xgblss = XGBoostLSS(
    SplineFlow(target_support="real",     # Specifies the support of the target. Options are "real", "positive", "positive_integer" or "unit_interval"
               count_bins=8,              # The number of segments comprising the spline.
               bound=bound,               # By adjusting the value, you can control the size of the bounding box and consequently control the range of inputs that the spline transform operates on.
               order="linear",            # The order of the spline. Options are "linear" or "quadratic".
               stabilization="None",      # Options are "None", "MAD" or "L2".
               loss_fn="nll"              # Loss function. Options are "nll" (negative log-likelihood) or "crps"(continuous ranked probability score).
            )    
)
```

--------------------------------

### Import required modules for multivariate distribution selection

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/How_To_Select_A_Multivariate_Distribution.ipynb

Load necessary distribution classes and utility functions for multivariate analysis.

```python
from xgboostlss.distributions.MVN import * 
from xgboostlss.distributions.MVT import *
from xgboostlss.distributions.MVN_LoRa import *

from xgboostlss.distributions.multivariate_distribution_utils import Multivariate_DistributionClass
from xgboostlss.datasets.data_loader import load_simulated_multivariate_gaussian_data
from sklearn import datasets
from sklearn.model_selection import train_test_split
```

--------------------------------

### Perform Hyper-Parameter Optimization

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/MVN_LowRank.ipynb

Initiate hyper-parameter optimization using `xgboostlss.hyper_opt`. Configure parameters like number of boosting rounds, cross-validation folds, early stopping, time budget, number of trials, and verbosity.

```python
np.random.seed(123)\nopt_param = xgblss.hyper_opt(param_dict,\n                             dtrain,\n                             num_boost_round=100,        # Number of boosting iterations.\n                             nfold=5,                    # Number of cv-folds.\n                             early_stopping_rounds=20,   # Number of early-stopping rounds\n                             max_minutes=120,            # Time budget in minutes, i.e., stop study after the given number of minutes.\n                             n_trials=20,                # The number of trials. If this argument is set to None, there is no limitation on the number of trials.\n                             silence=False,              # Controls the verbosity of the trail, i.e., user can silence the outputs of the trail.\n                             seed=123,                   # Seed used to generate cv-folds.\n                             hp_seed=None                # Seed for random number generator used in the Bayesian hyperparameter search.\n                            )
```

--------------------------------

### Initialize XGBoostLSS with Gaussian Mixture

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/GaussianMixture_Regression_CaliforniaHousing.ipynb

Configures the model using a Gaussian mixture distribution with specified components and hessian calculation mode.

```python
xgblss = XGBoostLSS(
    Mixture(
        Gaussian(response_fn="softplus"), 
        M = 4,
        tau=1.0,
        hessian_mode="individual",
    )
)
```

--------------------------------

### Preparing data for XGBoostLSS

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Gaussian_Regression.ipynb

Loads simulated Gaussian data and converts it into DMatrix format for model training.

```python
train, test = load_simulated_gaussian_data()
n_cpu = multiprocessing.cpu_count()

X_train, y_train = train.filter(regex="x"), train["y"].values
X_test, y_test = test.filter(regex="x"), test["y"].values

dtrain = xgb.DMatrix(X_train, label=y_train, nthread=n_cpu)
dtest = xgb.DMatrix(X_test, nthread=n_cpu)
```

--------------------------------

### Prepare simulated multivariate data

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/How_To_Select_A_Multivariate_Distribution.ipynb

Load and split the simulated dataset into training and testing sets.

```python
# Sample only 10% for plotting reasons
data = load_simulated_multivariate_gaussian_data().sample(frac=0.1, random_state=123).reset_index(drop=True).copy()

X, y = data["x"], data.filter(regex="y").values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=123)
```

--------------------------------

### Calculate and Prepare Quantiles for Plotting

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Gaussian_Regression.ipynb

Computes actual quantiles based on test data and prepares predicted quantiles for visualization.

```python
np.random.seed(123)

###
# Actual Quantiles
###
q1 = norm.ppf(quant_sel[0], loc = 10, scale = 1 + 4*((0.3 < test["x_true"].values) & (test["x_true"].values < 0.5)) + 2*(test["x_true"].values > 0.7))
q2 = norm.ppf(quant_sel[1], loc = 10, scale = 1 + 4*((0.3 < test["x_true"].values) & (test["x_true"].values < 0.5)) + 2*(test["x_true"].values > 0.7))
test["quant"] = np.where(test["y"].values < q1, 0, np.where(test["y"].values < q2, 1, 2))
test["alpha"] = np.where(test["y"].values <= q1, 1, np.where(test["y"].values >= q2, 1, 0))
df_quantiles = test[test["alpha"] == 1]

# Lower Bound
yl = list(set(q1))
yl.sort()
yl = [yl[2],yl[0],yl[2],yl[1],yl[1]]
sfunl = pd.DataFrame({"x_true":[0, 0.3, 0.5, 0.7, 1], "y":yl})

# Upper Bound
yu = list(set(q2))
yu.sort()
yu = [yu[0],yu[2],yu[0],yu[1],yu[1]]
sfunu = pd.DataFrame({"x_true":[0, 0.3, 0.5, 0.7, 1], "y":yu})

###
# Predicted Quantiles
###
test["lb"] = pred_quantiles.iloc[:,0]
test["ub"] = pred_quantiles.iloc[:,1]

###
# Plot
###
```

--------------------------------

### Load and prepare data for Dirichlet distribution

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Dirichlet_Regression.ipynb

Loads the Articlake dataset and prepares it for training. The features are separated from the target variables, and an xgb.DMatrix is created. This is suitable for compositional data where targets sum to a constant.

```python
data = load_articlake_data()

# Due to the small sample size, we train on the entire dataset
x_train = data[["depth"]]
y_train = data.drop(columns="depth").values
n_targets = y_train.shape[1]
dtrain = xgb.DMatrix(x_train, label=y_train, nthread=n_cpu)
```

--------------------------------

### Create and Plot Gaussian-Mixture Component Distributions

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/GaussianMixture_Regression_CaliforniaHousing.ipynb

Creates a mixture distribution and plots the component distributions. Requires `xgboostlss`, `torch`, `pandas`, and `seaborn`.

```python
torch.manual_seed(123)
mix_dist = xgblss.dist.create_mixture_distribution(mix_params)
gauss_dist = mix_dist._component_distribution
gauss_samples = pd.DataFrame(
    gauss_dist.sample((y_test.shape[0],)).reshape(-1, xgblss.dist.M).numpy(),
    columns = [f"Density {i+1}" for i in range(xgblss.dist.M)]
)

plt.figure(figsize=figure_size)
sns.kdeplot(gauss_samples, lw=2.5)
plt.title("Gaussian-Mixture Component Distributions", fontsize=20)
plt.show()
```

--------------------------------

### Distribution Selection with dist_select()

Source: https://context7.com/statmixedml/xgboostlss/llms.txt

Compares candidate distributions to identify the best fit based on negative log-likelihood.

```python
from xgboostlss.distributions.Gaussian import Gaussian
from xgboostlss.distributions.StudentT import StudentT
from xgboostlss.distributions.Gamma import Gamma
from xgboostlss.distributions.LogNormal import LogNormal
from xgboostlss.distributions import Beta, Cauchy, Weibull, Gumbel, Laplace
import numpy as np

# Generate sample data
np.random.seed(123)
y = np.random.gamma(shape=2, scale=2, size=1000)

# Define candidate distributions
candidate_distributions = [
    Gaussian,
    StudentT,
    Gamma,
    LogNormal,
    Weibull
]

# Use any distribution instance to call dist_select
dist = Gaussian()

# Compare distributions (set plot=True to visualize best fit)
fit_results = dist.dist_select(
    target=y,
    candidate_distributions=candidate_distributions,
    max_iter=100,
    plot=True,          # Requires matplotlib and seaborn
    figure_size=(10, 5)
)

print(fit_results)
# Output ranked by NLL (lower is better):
#       nll distribution
# 1  2834.12        Gamma
# 2  2901.45     LogNormal
# 3  2989.23      Gaussian
```

--------------------------------

### Import necessary libraries

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/MVN_Cholesky.ipynb

Imports all required libraries for data loading, manipulation, and model training.

```python
from xgboostlss.model import *
from xgboostlss.distributions.MVN import *
from xgboostlss.datasets.data_loader import load_simulated_multivariate_gaussian_data

from sklearn.model_selection import train_test_split
import pandas as pd
import multiprocessing
import plotnine
from plotnine import *
plotnine.options.figure_size = (15, 8)
n_cpu = multiprocessing.cpu_count()
```

--------------------------------

### Import necessary libraries

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Dirichlet_Regression.ipynb

Imports all required libraries for XGBoostLSS, data manipulation, and plotting.

```python
from xgboostlss.model import *
from xgboostlss.distributions.Dirichlet import *
from xgboostlss.datasets.data_loader import load_articlake_data

from sklearn.model_selection import train_test_split
import pandas as pd
import multiprocessing
import plotnine
from plotnine import *
plotnine.options.figure_size = (18, 9)
n_cpu = multiprocessing.cpu_count()
```

--------------------------------

### Visualize Distributional Parameters

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Gaussian_Regression.ipynb

Compares true Gaussian parameters against those predicted by the XGBoostLSS model.

```python
pred_params["x_true"] = X_test["x_true"].values
dist_params = list(xgblss.dist.param_dict.keys())

# Data with actual values
plot_df_actual = pd.melt(test[["x_true"] + dist_params],
                         id_vars="x_true",
                         value_vars=dist_params)
plot_df_actual["type"] = "TRUE"

# Data with predicted values
plot_df_predt = pd.melt(pred_params[["x_true"] + dist_params],
                        id_vars="x_true",
                        value_vars=dist_params)
plot_df_predt["type"] = "PREDICT"

plot_df = pd.concat([plot_df_predt, plot_df_actual])

plot_df["variable"] = plot_df.variable.str.upper()
plot_df["type"] = pd.Categorical(plot_df["type"], categories = ["PREDICT", "TRUE"])

(ggplot(plot_df,
        aes(x="x_true",
            y="value",
            color="type")) +
 geom_line(size=1.1) + 
 facet_wrap("variable",
            scales="free") + 
 labs(title="Parameters of univariate Gaussian predicted with XGBoostLSS",
      x="",
      y="") + 
 theme_bw(base_size=15) + 
 theme(legend_position="bottom",
       plot_title = element_text(hjust = 0.5),
       legend_title = element_blank())
)
```

--------------------------------

### Configure ZAGamma Distribution

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/ZAGamma_Regression.ipynb

Instantiates the XGBoostLSS model with the Zero-Adjusted Gamma distribution and specified loss and response functions.

```python
# Specifies Zero-Adjusted Gamma distribution. See ?ZAGamma for an overview.
xgblss = XGBoostLSS(
    ZAGamma(stabilization="None",        # Options are "None", "MAD", "L2".
            response_fn="exp",           # Function to transform the concentration and rate parameters, e.g., "exp" or "softplus".
            loss_fn="nll"                # Loss function. Options are "nll" (negative log-likelihood) or "crps"(continuous ranked probability score).)      
           )           
)
```

--------------------------------

### Configure Probabilistic Prediction Settings

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Dirichlet_Regression.ipynb

Sets the random seed for reproducibility and defines parameters for drawing samples and calculating quantiles from the predicted distribution.

```python
# Set seed for reproducibility
torch.manual_seed(123)

# Number of samples to draw from predicted distribution
n_samples = 1000
quant_sel = [0.05, 0.95] # Quantiles to calculate from predicted distribution
```

--------------------------------

### Prepare Simulation Data

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/ZAGamma_Regression.ipynb

Generates synthetic data and splits it into training and testing sets for XGBoost DMatrix format.

```python
# The simulation example closely follows https://towardsdatascience.com/zero-inflated-regression-c7dfc656d8af
np.random.seed(123)
n_samples = 1000

data = pd.DataFrame({"age": np.random.randint(1, 100, size=n_samples)})
data["income"] = np.where((data.age > 17) & (data.age < 70), 1500*data.age + 5000 + 10000*np.random.randn(n_samples), 0) / 1000

y = data["income"].values.reshape(-1,1) 
X = data.drop(columns="income")
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)

dtrain = xgb.DMatrix(X_train, label=y_train, nthread=n_cpu)
dtest = xgb.DMatrix(X_test, nthread=n_cpu)
```

--------------------------------

### Initialize XGBoostLSS with Gamma Distribution

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Gamma_Regression_CaliforniaHousing.ipynb

Configures the XGBoostLSS model to use the Gamma distribution. Specifies the response function as 'exp' and enables L2 stabilization for gradients and Hessians. The loss function is set to negative log-likelihood (nll).

```python
# Specifies Gamma distribution with exp response function and option to stabilize Gradient/Hessian. Type ?Gamma for an overview.
xgblss = XGBoostLSS(
    Gamma(stabilization="L2",     # Options are "None", "MAD", "L2".
          response_fn="exp",      # Function to transform the concentration and rate parameters, e.g., "exp" or "softplus".
          loss_fn="nll"           # Loss function. Options are "nll" (negative log-likelihood) or "crps"(continuous ranked probability score).
         )
)
```

--------------------------------

### Univariate Distributions

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/distributions.md

Details on univariate distributions, including their PyTorch distribution links, XGBoostlss class names, types, output ranges, and parameter counts.

```APIDOC
## Univariate Distributions

### Mixture
- **Description**: Continuous & Discrete Count (Univariate)
- **Output Range**: $y \in (-\infty,\infty)$, $y \in [0, \infty)$, $y \in [0, 1]$, $y \in (0, 1, 2, 3, \ldots)$
- **Parameters**: CompDist + M

### Negative Binomial
- **Description**: Discrete Count (Univariate)
- **Output Range**: $y \in (0, 1, 2, 3, \ldots)$
- **Parameters**: 2

### Poisson
- **Description**: Discrete Count (Univariate)
- **Output Range**: $y \in (0, 1, 2, 3, \ldots)$
- **Parameters**: 1

### Spline Flow
- **Description**: Continuous & Discrete Count (Univariate)
- **Output Range**: $y \in (-\infty,\infty)$, $y \in [0, \infty)$, $y \in [0, 1]$, $y \in (0, 1, 2, 3, \ldots)$
- **Parameters**: 2xcount_bins + (count_bins-1) (order=quadratic) or 3xcount_bins + (count_bins-1) (order=linear)

### Student-T
- **Description**: Continuous (Univariate)
- **Output Range**: $y \in (-\infty,\infty)$
- **Parameters**: 3

### Weibull
- **Description**: Continuous (Univariate)
- **Output Range**: $y \in [0, \infty)$
- **Parameters**: 2

### Zero-Adjusted Beta
- **Description**: Discrete-Continuous (Univariate)
- **Output Range**: $y \in [0, 1)$
- **Parameters**: 3

### Zero-Adjusted Gamma
- **Description**: Discrete-Continuous (Univariate)
- **Output Range**: $y \in [0, \infty)$
- **Parameters**: 3

### Zero-Adjusted LogNormal
- **Description**: Discrete-Continuous (Univariate)
- **Output Range**: $y \in [0, \infty)$
- **Parameters**: 3

### Zero-Inflated Negative Binomial
- **Description**: Discrete-Count (Univariate)
- **Output Range**: $y \in [0, 1, 2, 3, \ldots)$
- **Parameters**: 3
```

--------------------------------

### Define and Select Candidate SplineFlow Models

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/SplineFlow_Regression.ipynb

Configures a list of SplineFlow instances with different bin counts and orders, then executes a selection process to find the optimal flow for the training data.

```python
bound = np.max([np.abs(y_train.min()), y_train.max()])
target_support = "real"

candidate_flows = [

    SplineFlow(target_support=target_support, count_bins=2,  bound=bound, order="linear"),
    SplineFlow(target_support=target_support, count_bins=4,  bound=bound, order="linear"),
    SplineFlow(target_support=target_support, count_bins=6,  bound=bound, order="linear"),
    SplineFlow(target_support=target_support, count_bins=8,  bound=bound, order="linear"),
    SplineFlow(target_support=target_support, count_bins=12, bound=bound, order="linear"),
    SplineFlow(target_support=target_support, count_bins=16, bound=bound, order="linear"),
    SplineFlow(target_support=target_support, count_bins=20, bound=bound, order="linear"),

    SplineFlow(target_support=target_support, count_bins=2,  bound=bound, order="quadratic"),
    SplineFlow(target_support=target_support, count_bins=4,  bound=bound, order="quadratic"),
    SplineFlow(target_support=target_support, count_bins=6,  bound=bound, order="quadratic"),
    SplineFlow(target_support=target_support, count_bins=8,  bound=bound, order="quadratic"),
    SplineFlow(target_support=target_support, count_bins=12, bound=bound, order="quadratic"),
    SplineFlow(target_support=target_support, count_bins=16, bound=bound, order="quadratic"),
    SplineFlow(target_support=target_support, count_bins=20, bound=bound, order="quadratic"),
    
] 

flow_nll = NormalizingFlowClass().flow_select(target=y_train, candidate_flows=candidate_flows, max_iter=50, plot=True, figure_size=(12, 5))
flow_nll
```

--------------------------------

### Import necessary libraries

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/How_To_Select_A_Univariate_Distribution.ipynb

Imports required modules for distributions, data handling, and model splitting.

```python
from xgboostlss.distributions import *
from xgboostlss.distributions.distribution_utils import DistributionClass
from sklearn import datasets
from sklearn.model_selection import train_test_split
```

--------------------------------

### Inspect Prediction Samples

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/MVT_Cholesky.ipynb

View the head and tail of the generated prediction samples.

```python
pred_samples.head()
```

```python
pred_samples.tail()
```

--------------------------------

### Configure XGBoostLSS for Multivariate Student-T Regression

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/MVT_Cholesky.ipynb

Initializes the XGBoostLSS model with a Multivariate Student-T distribution. It specifies the number of targets (D), the stabilization method for the covariance matrix, the response function for the Cholesky factor transformation, and the loss function.

```python
# Specifies a multivariate Student-T distribution, using the Cholesky decompoisition. See ?MVT for details.
xgblss = XGBoostLSS(
    MVT(D=n_targets,                      # Specifies the number of targets
        stabilization="None",             # Options are "None", "MAD", "L2".
        response_fn="exp",                # Function to transform the lower-triangular factor of the covariance, e.g., "exp" or "softplus".
        loss_fn="nll"                     # Loss function, i.e., nll.
       )               
)

```

--------------------------------

### Train XGBoostLSS Model with train()

Source: https://context7.com/statmixedml/xgboostlss/llms.txt

Trains an XGBoostLSS model using standard XGBoost parameters. The framework automatically handles multi-parameter optimization. Ensure data is loaded and prepared into DMatrix format.

```python
from xgboostlss.model import XGBoostLSS
from xgboostlss.distributions.Gaussian import Gaussian
from xgboostlss.datasets.data_loader import load_simulated_gaussian_data
import xgboost as xgb

# Load example data
train, test = load_simulated_gaussian_data()
X_train, y_train = train.filter(regex="x"), train["y"].values
X_test, y_test = test.filter(regex="x"), test["y"].values

dtrain = xgb.DMatrix(X_train, label=y_train)
dtest = xgb.DMatrix(X_test)

# Initialize model
xgblss = XGBoostLSS(Gaussian(stabilization="None", response_fn="exp", loss_fn="nll"))

# Training parameters
params = {
    "eta": 0.1,                # Learning rate
    "max_depth": 3,            # Maximum tree depth
    "min_child_weight": 1,     # Minimum sum of instance weight in child
    "subsample": 0.9,          # Subsample ratio
    "colsample_bytree": 0.9,   # Feature subsample ratio
    "booster": "gbtree"        # Booster type
}

# Train the model
xgblss.train(
    params=params,
    dtrain=dtrain,
    num_boost_round=100,
    evals=[(dtrain, "train"), (dtest, "test")],
    early_stopping_rounds=10,
    verbose_eval=10
)
```

--------------------------------

### Initialize XGBoostLSS with Dirichlet distribution

Source: https://github.com/statmixedml/xgboostlss/blob/master/docs/examples/Dirichlet_Regression.ipynb

Initializes the XGBoostLSS model with the Dirichlet distribution. Specify the number of targets (D), stabilization method, response function, and loss function. The Dirichlet distribution is used for compositional data.

```python
# Specifies a Dirichlet distribution. See ?Dirichlet for details.
xgblss = XGBoostLSS(
    Dirichlet(D=n_targets,                      # Specifies the number of targets
              stabilization="None",             # Options are "None", "MAD", "L2".
              response_fn="exp",                # Function to transform the lower-triangular factor of the covariance, e.g., "exp", "relu" or "softplus".
              loss_fn="nll"                     # Loss function, i.e., nll.
             )               
)
```