求助：构建依赖Laplace变量的无NaN截断Laplace/正态分布

Problem Breakdown #

Your current workflow generates full-range Laplace noise, adds it to your linear prediction, then discards out-of-bounds values by setting them to NaNs. This leads to data loss and defeats the goal of having every x map to a valid, bounded y. The fix is to generate noise that strictly stays within your constraints from the start.

Solution: Direct Truncated Laplace Sampling #

Instead of post-processing with NaNs, we can generate truncated Laplace noise tailored to each x's bounds. This uses the inverse CDF method to efficiently produce samples that exactly follow the truncated Laplace distribution, ensuring no invalid values are created in the first place.

Full Modified Code #

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import gaussian_kde, truncnorm

slope = 0.2237
intercept = 1.066
spread = 4.8719

def truncated_laplace(loc, scale, lower, upper, size):
    """Generate samples from a truncated Laplace distribution."""
    def laplace_cdf(z, loc, scale):
        """CDF for the Laplace distribution (loc = center, scale = spread)."""
        return 0.5 * (1 + np.sign(z - loc) * (1 - np.exp(-np.abs(z - loc)/scale)))
    
    # Calculate CDF values at our truncation bounds
    cdf_lower = laplace_cdf(lower, loc, scale)
    cdf_upper = laplace_cdf(upper, loc, scale)
    
    # Generate uniform samples between the two CDF values
    u = np.random.uniform(cdf_lower, cdf_upper, size=size)
    
    # Invert the CDF to get truncated Laplace samples
    z = np.empty(size)
    # Handle right half of the Laplace distribution (samples >= loc)
    mask_right = u >= 0.5
    z[mask_right] = loc - scale * np.log(2 * (1 - u[mask_right]))
    # Handle left half of the Laplace distribution (samples < loc)
    mask_left = ~mask_right
    z[mask_left] = loc + scale * np.log(2 * u[mask_left])
    
    return z

def dependency(x):
    y_lin = slope * x + intercept
    # Compute noise bounds relative to the linear prediction y_lin
    noise_lower = (slope / spread * 3 * x) - y_lin
    noise_upper = (slope * spread / 3 * x + 2 * intercept) - y_lin
    
    # Generate truncated noise tailored to each x's bounds
    noise = truncated_laplace(loc=0, scale=spread, lower=noise_lower, upper=noise_upper, size=len(x))
    
    # Calculate final y values
    y = y_lin + noise
    
    # Optional validation to ensure no values are out of bounds
    assert np.all(y >= slope / spread * 3 * x), "y values below lower constraint!"
    assert np.all(y <= slope * spread / 3 * x + 2 * intercept), "y values above upper constraint!"
    return y

# Generate x from truncated normal distribution
max_x = 100
min_x = 1
mean_x = 40
sigma_x = 25
x = truncnorm((min_x - mean_x)/sigma_x, (max_x - mean_x)/sigma_x, loc=mean_x, scale=sigma_x).rvs(5000)

# Generate y with no NaNs
y = dependency(x)

# Plotting (no need to filter NaNs anymore!)
xx = np.linspace(min_x, max_x, 100)
yy = slope * xx + intercept
lower = slope/spread*3*xx
upper = slope*spread/3*xx + 2*intercept

xy = np.vstack([x, y])
z = gaussian_kde(xy)(xy)
idz = z.argsort()
x_sorted, y_sorted, z_sorted = x[idz], y[idz], z[idz]

fig, ax = plt.subplots(figsize=(5, 5))
plt.plot(xx, upper, 'r-.', label='upper constraint')
plt.plot(xx, lower, 'r--', label='lower constraint')
ax.scatter(x_sorted, y_sorted, c=z_sorted, s=3)
plt.xlabel(r'$\bf X_{truncated\ normal}$')
plt.ylabel(r'$\bf Y_{truncated\ Laplace}$')
plt.plot(xx, yy, 'r', label='regression model')
plt.legend()
plt.tight_layout()
plt.show()

Key Improvements #

• No NaNs: Every x gets a valid y within your specified bounds—no data filtering required.
• True Truncated Laplace: The sampling method ensures noise follows the exact truncated Laplace distribution, not just a filtered version of the full distribution.
• Validation Checks: The assert statements confirm your constraints are respected (you can remove them in production for speed).

Quick Note on Bounds #

We verified your constraints are valid for all x in [1, 100]: given your spread value (~4.87), lower(x) < y_lin < upper(x) always holds, so the truncation range is meaningful for every input.

内容的提问来源于stack exchange，提问作者M. Weeker