distribution_fitting

distribution_fitting

Module: distribution_fitting.py

This module contains functions for statistical distribution fitting and quantile mapping used in the bias correction workflow. It includes functions for: - Calculating L-moments and L-moment ratios (using unbiased PWM estimators from Hosking 1990). - Fitting a gamma distribution using MLE (scipy.stats.gamma.fit with floc=0). - Fitting a Generalized Pareto Distribution (GPD) to data above a threshold. - Performing K-Fold cross-validation to obtain stable GPD parameters. - Applying gamma distribution-based quantile mapping with GPD tail adjustment using proper conditional probability mapping for extreme values.

Dependencies: - Imports default parameters from the config module. - Uses scipy.stats and sklearn.model_selection for statistical calculations.

Author: Benny Istanto Applied Climatology Study Program, Department of Geophysics and Meteorology, Bogor Agricultural University, Indonesia Email: bennyistanto@apps.ipb.ac.id

with supervision from Prof. Rizaldi Boer and Dr. I Putu Santikayasa

Update: 2026.03

Functions

Name	Description
calculate_l_moments	Calculate L-moments and L-moment ratios for the given data.
cross_validate_gpd	Cross-validate GPD fitting by splitting data into folds.
fit_cpc_parameters_on_native_grid	Fit gamma and GPD distribution parameters at each CPC native-resolution cell.
fit_gamma_distribution	Fit a gamma distribution to the data using Maximum Likelihood Estimation (MLE)
fit_generalized_pareto_distribution	Fit a Generalized Pareto Distribution (GPD) to the excesses above the threshold.
gamma_quantile_mapping	Apply gamma distribution-based quantile mapping with tail adjustment
gamma_quantile_mapping_precomputed	Apply gamma distribution-based quantile mapping with GPD tail adjustment,
interpolate_cpc_params_to_imerg_grid	Bilinearly interpolate CPC distribution parameters from native ~0.5° to

calculate_l_moments

distribution_fitting.calculate_l_moments(data)

Calculate L-moments and L-moment ratios for the given data.

L-moments are statistics used to describe the shape of a probability distribution. They are analogous to conventional moments but can be more robust to outliers.

Uses the unbiased PWM estimators from Hosking (1990): b_r = (1/n) * sum_{i=1}^{n} C(i-1,r)/C(n-1,r) * x_{i:n} where x_{i:n} are the order statistics.

Parameters: data (numpy.ndarray): The input data.

Returns: tuple: L-moments (l1, l2, l3, l4) and L-moment ratios (t2, t3, t4)

cross_validate_gpd

distribution_fitting.cross_validate_gpd(
    data,
    threshold,
    n_splits=N_SPLITS_GPD_CROSSVALIDATE,
)

Cross-validate GPD fitting by splitting data into folds.

This function uses K-Fold cross-validation to assess the stability and reliability of the GPD parameter estimates.

Parameters: data (numpy.ndarray): Array of data values. threshold (float): Threshold value for defining the excesses. n_splits (int, optional): Number of cross-validation splits. Default is 5.

Returns: tuple: Averaged parameters of the GPD from cross-validation (shape, location, scale).

fit_cpc_parameters_on_native_grid

distribution_fitting.fit_cpc_parameters_on_native_grid(cpc_native_dekad)

Fit gamma and GPD distribution parameters at each CPC native-resolution cell.

This avoids redundantly fitting the same CPC time series 25 times (as happens when CPC is nearest-neighbour regridded to 0.1°). Instead, parameters are fitted once per ~0.5° cell and later interpolated to the IMERG grid.

Based on the BCSD principle (Wood et al. 2004): correct at reference resolution, then disaggregate smoothly.

Parameters

Name	Type	Description	Default
cpc_native_dekad	xarray.DataArray	CPC precipitation at native ~0.5° resolution for one dekad across all years. Shape (n_time, n_lat_cpc, n_lon_cpc).	required

Returns

Name	Type	Description
	dict of xarray.DataArray	Dictionary with 9 parameter arrays at CPC native resolution: ‘gamma_shape’, ‘gamma_scale’, ‘gpd_threshold’, ‘gpd_shape’, ‘gpd_loc’, ‘gpd_scale’, ‘upper_cap’, ‘p_threshold’, ‘p_dry_cpc’. Ocean/invalid cells are NaN.

Notes

Fixes applied (2026.04): * Fix C (Cannon 2015 §3.2): gamma distribution fitted on WET-day values only (> WET_DAY_THRESHOLD), not on all values including zeros. The GPD threshold is also computed from wet values only. * Fix A (Coles 2001): GPD fitting inside cross_validate_gpd uses floc=0 (see fit_generalized_pareto_distribution). * New: p_dry_cpc is returned per cell so the downstream mapping step can apply Cannon dry-day handling using the interpolated dry-day frequency at each IMERG pixel.

fit_gamma_distribution

distribution_fitting.fit_gamma_distribution(data)

Fit a gamma distribution to the data using Maximum Likelihood Estimation (MLE) via scipy.stats.gamma.fit with location fixed at zero.

For precipitation data, the gamma distribution is a natural choice as it models positive-valued, right-skewed random variables. Fixing loc=0 ensures the distribution starts at zero, which is physically meaningful for precipitation.

Parameters: data (numpy.ndarray): Array of data values.

Returns: tuple: Fitted parameters (shape, loc, scale) of the gamma distribution.

fit_generalized_pareto_distribution

distribution_fitting.fit_generalized_pareto_distribution(data, threshold)

Fit a Generalized Pareto Distribution (GPD) to the excesses above the threshold.

The GPD is often used in extreme value theory to model the tail of a distribution. It’s particularly useful for modeling events that exceed a high threshold.

Parameters: data (numpy.ndarray): Array of data values. threshold (float): Threshold value for defining the excesses.

Returns: tuple: Fitted parameters of the GPD (shape, location, scale).

gamma_quantile_mapping

distribution_fitting.gamma_quantile_mapping(imerg_values, cpc_values)

Apply gamma distribution-based quantile mapping with tail adjustment to correct the distribution of precipitation data.

This function fits gamma distributions to the IMERG and CPC precipitation values using MLE (scipy.stats.gamma.fit). It then computes the cumulative distribution function (CDF) of the IMERG values and applies the inverse CDF of the CPC values to obtain the corrected precipitation values. Additionally, it adjusts the tails using the Generalized Pareto Distribution (GPD) to better capture extreme values.

The GPD tail adjustment uses proper conditional probability mapping: For values above the threshold, the unconditional CDF value is mapped to the conditional exceedance probability before applying the GPD inverse CDF.

Parameters: imerg_values (numpy.ndarray): Array of IMERG precipitation values. cpc_values (numpy.ndarray): Array of CPC precipitation values.

Returns: numpy.ndarray: Corrected precipitation values after gamma quantile mapping with tail adjustment.

Notes

Fixes applied (2026.04): * Fix C (Cannon 2015 §3.2): gamma distributions are fitted on wet-day values only (> WET_DAY_THRESHOLD). The dry-day fraction is handled explicitly via the unconditional CDF construction, and sub-threshold IMERG values are mapped either to zero (when the unconditional probability is below the CPC dry-day fraction) or through the conditional wet CDF of CPC. This prevents drizzle days from contaminating the gamma body and matches the CPC dry-day frequency. * Fix B: GPD substitution is triggered on IMERG’s own 80th percentile of wet values, not on the CPC threshold, and the conditional probability is computed consistently in unconditional CDF space. * Fix A (Coles 2001): GPD fitting uses floc=0 internally.

gamma_quantile_mapping_precomputed

distribution_fitting.gamma_quantile_mapping_precomputed(
    imerg_values,
    cpc_gamma_shape,
    cpc_gamma_scale,
    cpc_gpd_threshold,
    cpc_gpd_shape,
    cpc_gpd_loc,
    cpc_gpd_scale,
    cpc_upper_cap,
    cpc_p_threshold,
    p_dry_cpc,
)

Apply gamma distribution-based quantile mapping with GPD tail adjustment, using pre-computed CPC-side parameters instead of fitting CPC per pixel.

This is the BCSD-path counterpart to gamma_quantile_mapping(). Statistical logic is identical, but CPC parameters come from the native-resolution fit (fit_cpc_parameters_on_native_grid) and are bilinearly interpolated to the IMERG grid, eliminating 0.5° block artefacts.

Parameters

Name	Type	Description	Default
imerg_values	numpy.ndarray	IMERG precipitation time series for one pixel.	required
cpc_gamma_shape	float	Pre-computed CPC gamma shape (fitted on wet values only, Fix C).	required
cpc_gamma_scale	float	Pre-computed CPC gamma scale (fitted on wet values only, Fix C).	required
cpc_gpd_threshold	float	Pre-computed CPC GPD threshold (80th percentile of wet values).	required
cpc_gpd_shape	float	Pre-computed CPC GPD shape (fitted with floc=0, Fix A).	required
cpc_gpd_loc	float	Pre-computed CPC GPD location.	required
cpc_gpd_scale	float	Pre-computed CPC GPD scale.	required
cpc_upper_cap	float	Pre-computed CPC upper cap (99.9th percentile of wet values).	required
cpc_p_threshold	float	Pre-computed CPC gamma CDF value at the GPD threshold.	required
p_dry_cpc	float	Pre-computed CPC dry-day fraction (fraction of days <= WET_DAY_THRESHOLD). Required for Cannon 2015 §3.2 dry-day handling (Fix C).	required

Returns

Name	Type	Description
	numpy.ndarray	Corrected precipitation values after quantile mapping.

Notes

Fixes applied (2026.04): * Fix C (Cannon 2015 §3.2): IMERG gamma fitted on wet-day values only; Cannon dry-day handling uses p_dry_imerg (per pixel) and p_dry_cpc (interpolated from native CPC fit) to match dry-day frequencies. * Fix B: GPD substitution triggered on IMERG-side 80th percentile of wet values, with conditional probability computed in the unconditional CDF space for consistency between gamma body and GPD tail. * Fix A (Coles 2001): inherited via cpc_gpd_* which were fitted with floc=0 upstream in cross_validate_gpd.

interpolate_cpc_params_to_imerg_grid

distribution_fitting.interpolate_cpc_params_to_imerg_grid(
    cpc_params,
    target_lat,
    target_lon,
)

Bilinearly interpolate CPC distribution parameters from native ~0.5° to the IMERG 0.1° grid. Uses a two-pass approach: bilinear first, then nearest-neighbour to fill boundary NaN values.

Parameters

Name	Type	Description	Default
cpc_params	dict of xarray.DataArray	CPC distribution parameters at native resolution, as returned by fit_cpc_parameters_on_native_grid().	required
target_lat	numpy.ndarray	Target latitude values (IMERG grid).	required
target_lon	numpy.ndarray	Target longitude values (IMERG grid).	required

Returns

Name	Type	Description
	dict of xarray.DataArray	Interpolated parameters at IMERG resolution.

# distribution_fitting { #src.distribution_fitting } `distribution_fitting` Module: distribution_fitting.py This module contains functions for statistical distribution fitting and quantile mapping used in the bias correction workflow. It includes functions for: - Calculating L-moments and L-moment ratios (using unbiased PWM estimators from Hosking 1990). - Fitting a gamma distribution using MLE (scipy.stats.gamma.fit with floc=0). - Fitting a Generalized Pareto Distribution (GPD) to data above a threshold. - Performing K-Fold cross-validation to obtain stable GPD parameters. - Applying gamma distribution-based quantile mapping with GPD tail adjustment using proper conditional probability mapping for extreme values. Dependencies: - Imports default parameters from the config module. - Uses scipy.stats and sklearn.model_selection for statistical calculations. **Author**: Benny Istanto Applied Climatology Study Program, Department of Geophysics and Meteorology, Bogor Agricultural University, Indonesia Email: bennyistanto@apps.ipb.ac.id with supervision from Prof. Rizaldi Boer and Dr. I Putu Santikayasa Update: 2026.03 ## Functions | Name | Description | | --- | --- | | [calculate_l_moments](#src.distribution_fitting.calculate_l_moments) | Calculate L-moments and L-moment ratios for the given data. | | [cross_validate_gpd](#src.distribution_fitting.cross_validate_gpd) | Cross-validate GPD fitting by splitting data into folds. | | [fit_cpc_parameters_on_native_grid](#src.distribution_fitting.fit_cpc_parameters_on_native_grid) | Fit gamma and GPD distribution parameters at each CPC native-resolution cell. | | [fit_gamma_distribution](#src.distribution_fitting.fit_gamma_distribution) | Fit a gamma distribution to the data using Maximum Likelihood Estimation (MLE) | | [fit_generalized_pareto_distribution](#src.distribution_fitting.fit_generalized_pareto_distribution) | Fit a Generalized Pareto Distribution (GPD) to the excesses above the threshold. | | [gamma_quantile_mapping](#src.distribution_fitting.gamma_quantile_mapping) | Apply gamma distribution-based quantile mapping with tail adjustment | | [gamma_quantile_mapping_precomputed](#src.distribution_fitting.gamma_quantile_mapping_precomputed) | Apply gamma distribution-based quantile mapping with GPD tail adjustment, | | [interpolate_cpc_params_to_imerg_grid](#src.distribution_fitting.interpolate_cpc_params_to_imerg_grid) | Bilinearly interpolate CPC distribution parameters from native ~0.5° to | ### calculate_l_moments { #src.distribution_fitting.calculate_l_moments } ```python distribution_fitting.calculate_l_moments(data) ``` Calculate L-moments and L-moment ratios for the given data. L-moments are statistics used to describe the shape of a probability distribution. They are analogous to conventional moments but can be more robust to outliers. Uses the unbiased PWM estimators from Hosking (1990): b_r = (1/n) * sum_{i=1}^{n} C(i-1,r)/C(n-1,r) * x_{i:n} where x_{i:n} are the order statistics. Parameters: data (numpy.ndarray): The input data. Returns: tuple: L-moments (l1, l2, l3, l4) and L-moment ratios (t2, t3, t4) ### cross_validate_gpd { #src.distribution_fitting.cross_validate_gpd } ```python distribution_fitting.cross_validate_gpd( data, threshold, n_splits=N_SPLITS_GPD_CROSSVALIDATE, ) ``` Cross-validate GPD fitting by splitting data into folds. This function uses K-Fold cross-validation to assess the stability and reliability of the GPD parameter estimates. Parameters: data (numpy.ndarray): Array of data values. threshold (float): Threshold value for defining the excesses. n_splits (int, optional): Number of cross-validation splits. Default is 5. Returns: tuple: Averaged parameters of the GPD from cross-validation (shape, location, scale). ### fit_cpc_parameters_on_native_grid { #src.distribution_fitting.fit_cpc_parameters_on_native_grid } ```python distribution_fitting.fit_cpc_parameters_on_native_grid(cpc_native_dekad) ``` Fit gamma and GPD distribution parameters at each CPC native-resolution cell. This avoids redundantly fitting the same CPC time series 25 times (as happens when CPC is nearest-neighbour regridded to 0.1°). Instead, parameters are fitted once per ~0.5° cell and later interpolated to the IMERG grid. Based on the BCSD principle (Wood et al. 2004): correct at reference resolution, then disaggregate smoothly. #### Parameters {.doc-section .doc-section-parameters} | Name | Type | Description | Default | |------------------|----------------------------------------------------|--------------------------------------------------------------------------------------------------------------------|------------| | cpc_native_dekad | [xarray](`xarray`).[DataArray](`xarray.DataArray`) | CPC precipitation at native ~0.5° resolution for one dekad across all years. Shape (n_time, n_lat_cpc, n_lon_cpc). | _required_ | #### Returns {.doc-section .doc-section-returns} | Name | Type | Description | |--------|--------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------| | | dict of xarray.DataArray | Dictionary with 9 parameter arrays at CPC native resolution: 'gamma_shape', 'gamma_scale', 'gpd_threshold', 'gpd_shape', 'gpd_loc', 'gpd_scale', 'upper_cap', 'p_threshold', 'p_dry_cpc'. Ocean/invalid cells are NaN. | #### Notes {.doc-section .doc-section-notes} Fixes applied (2026.04): * Fix C (Cannon 2015 §3.2): gamma distribution fitted on WET-day values only (> WET_DAY_THRESHOLD), not on all values including zeros. The GPD threshold is also computed from wet values only. * Fix A (Coles 2001): GPD fitting inside cross_validate_gpd uses floc=0 (see fit_generalized_pareto_distribution). * New: p_dry_cpc is returned per cell so the downstream mapping step can apply Cannon dry-day handling using the interpolated dry-day frequency at each IMERG pixel. ### fit_gamma_distribution { #src.distribution_fitting.fit_gamma_distribution } ```python distribution_fitting.fit_gamma_distribution(data) ``` Fit a gamma distribution to the data using Maximum Likelihood Estimation (MLE) via scipy.stats.gamma.fit with location fixed at zero. For precipitation data, the gamma distribution is a natural choice as it models positive-valued, right-skewed random variables. Fixing loc=0 ensures the distribution starts at zero, which is physically meaningful for precipitation. Parameters: data (numpy.ndarray): Array of data values. Returns: tuple: Fitted parameters (shape, loc, scale) of the gamma distribution. ### fit_generalized_pareto_distribution { #src.distribution_fitting.fit_generalized_pareto_distribution } ```python distribution_fitting.fit_generalized_pareto_distribution(data, threshold) ``` Fit a Generalized Pareto Distribution (GPD) to the excesses above the threshold. The GPD is often used in extreme value theory to model the tail of a distribution. It's particularly useful for modeling events that exceed a high threshold. Parameters: data (numpy.ndarray): Array of data values. threshold (float): Threshold value for defining the excesses. Returns: tuple: Fitted parameters of the GPD (shape, location, scale). ### gamma_quantile_mapping { #src.distribution_fitting.gamma_quantile_mapping } ```python distribution_fitting.gamma_quantile_mapping(imerg_values, cpc_values) ``` Apply gamma distribution-based quantile mapping with tail adjustment to correct the distribution of precipitation data. This function fits gamma distributions to the IMERG and CPC precipitation values using MLE (scipy.stats.gamma.fit). It then computes the cumulative distribution function (CDF) of the IMERG values and applies the inverse CDF of the CPC values to obtain the corrected precipitation values. Additionally, it adjusts the tails using the Generalized Pareto Distribution (GPD) to better capture extreme values. The GPD tail adjustment uses proper conditional probability mapping: For values above the threshold, the unconditional CDF value is mapped to the conditional exceedance probability before applying the GPD inverse CDF. Parameters: imerg_values (numpy.ndarray): Array of IMERG precipitation values. cpc_values (numpy.ndarray): Array of CPC precipitation values. Returns: numpy.ndarray: Corrected precipitation values after gamma quantile mapping with tail adjustment. #### Notes {.doc-section .doc-section-notes} Fixes applied (2026.04): * Fix C (Cannon 2015 §3.2): gamma distributions are fitted on wet-day values only (> WET_DAY_THRESHOLD). The dry-day fraction is handled explicitly via the unconditional CDF construction, and sub-threshold IMERG values are mapped either to zero (when the unconditional probability is below the CPC dry-day fraction) or through the conditional wet CDF of CPC. This prevents drizzle days from contaminating the gamma body and matches the CPC dry-day frequency. * Fix B: GPD substitution is triggered on IMERG's own 80th percentile of wet values, not on the CPC threshold, and the conditional probability is computed consistently in unconditional CDF space. * Fix A (Coles 2001): GPD fitting uses floc=0 internally. ### gamma_quantile_mapping_precomputed { #src.distribution_fitting.gamma_quantile_mapping_precomputed } ```python distribution_fitting.gamma_quantile_mapping_precomputed( imerg_values, cpc_gamma_shape, cpc_gamma_scale, cpc_gpd_threshold, cpc_gpd_shape, cpc_gpd_loc, cpc_gpd_scale, cpc_upper_cap, cpc_p_threshold, p_dry_cpc, ) ``` Apply gamma distribution-based quantile mapping with GPD tail adjustment, using pre-computed CPC-side parameters instead of fitting CPC per pixel. This is the BCSD-path counterpart to gamma_quantile_mapping(). Statistical logic is identical, but CPC parameters come from the native-resolution fit (fit_cpc_parameters_on_native_grid) and are bilinearly interpolated to the IMERG grid, eliminating 0.5° block artefacts. #### Parameters {.doc-section .doc-section-parameters} | Name | Type | Description | Default | |-------------------|---------------------------------------------|------------------------------------------------------------------------------------------------------------------------------------|------------| | imerg_values | [numpy](`numpy`).[ndarray](`numpy.ndarray`) | IMERG precipitation time series for one pixel. | _required_ | | cpc_gamma_shape | [float](`float`) | Pre-computed CPC gamma shape (fitted on wet values only, Fix C). | _required_ | | cpc_gamma_scale | [float](`float`) | Pre-computed CPC gamma scale (fitted on wet values only, Fix C). | _required_ | | cpc_gpd_threshold | [float](`float`) | Pre-computed CPC GPD threshold (80th percentile of wet values). | _required_ | | cpc_gpd_shape | [float](`float`) | Pre-computed CPC GPD shape (fitted with floc=0, Fix A). | _required_ | | cpc_gpd_loc | [float](`float`) | Pre-computed CPC GPD location. | _required_ | | cpc_gpd_scale | [float](`float`) | Pre-computed CPC GPD scale. | _required_ | | cpc_upper_cap | [float](`float`) | Pre-computed CPC upper cap (99.9th percentile of wet values). | _required_ | | cpc_p_threshold | [float](`float`) | Pre-computed CPC gamma CDF value at the GPD threshold. | _required_ | | p_dry_cpc | [float](`float`) | Pre-computed CPC dry-day fraction (fraction of days <= WET_DAY_THRESHOLD). Required for Cannon 2015 §3.2 dry-day handling (Fix C). | _required_ | #### Returns {.doc-section .doc-section-returns} | Name | Type | Description | |--------|---------------------------------------------|--------------------------------------------------------| | | [numpy](`numpy`).[ndarray](`numpy.ndarray`) | Corrected precipitation values after quantile mapping. | #### Notes {.doc-section .doc-section-notes} Fixes applied (2026.04): * Fix C (Cannon 2015 §3.2): IMERG gamma fitted on wet-day values only; Cannon dry-day handling uses p_dry_imerg (per pixel) and p_dry_cpc (interpolated from native CPC fit) to match dry-day frequencies. * Fix B: GPD substitution triggered on IMERG-side 80th percentile of wet values, with conditional probability computed in the unconditional CDF space for consistency between gamma body and GPD tail. * Fix A (Coles 2001): inherited via cpc_gpd_* which were fitted with floc=0 upstream in cross_validate_gpd. ### interpolate_cpc_params_to_imerg_grid { #src.distribution_fitting.interpolate_cpc_params_to_imerg_grid } ```python distribution_fitting.interpolate_cpc_params_to_imerg_grid( cpc_params, target_lat, target_lon, ) ``` Bilinearly interpolate CPC distribution parameters from native ~0.5° to the IMERG 0.1° grid. Uses a two-pass approach: bilinear first, then nearest-neighbour to fill boundary NaN values. #### Parameters {.doc-section .doc-section-parameters} | Name | Type | Description | Default | |------------|---------------------------------------------|-------------------------------------------------------------------------------------------------------|------------| | cpc_params | dict of xarray.DataArray | CPC distribution parameters at native resolution, as returned by fit_cpc_parameters_on_native_grid(). | _required_ | | target_lat | [numpy](`numpy`).[ndarray](`numpy.ndarray`) | Target latitude values (IMERG grid). | _required_ | | target_lon | [numpy](`numpy`).[ndarray](`numpy.ndarray`) | Target longitude values (IMERG grid). | _required_ | #### Returns {.doc-section .doc-section-returns} | Name | Type | Description | |--------|--------------------------|----------------------------------------------| | | dict of xarray.DataArray | Interpolated parameters at IMERG resolution. |