实现Johnson SU分布的参数计算和优化过程
Johnson SU 分布参数计算与优化
实现Johnson SU分布的参数计算和优化过程。Johnson SU分布是一种灵活的分布族，可以拟合各种形状的数据分布。
1. 导入必要的库
import numpy as np
from scipy.optimize import minimize
from sklearn.linear_model import ElasticNet
from sklearn.model_selection import cross_val_score
from bayes_opt import BayesianOptimization
import warnings
warnings.filterwarnings('ignore')
# 设置随机种子保证结果可复现
np.random.seed(42)
2. 计算位置参数γ
def calculate_gamma(x):"""计算位置参数γ参数:x -- 输入数据数组返回:gamma -- 位置参数γ"""# ⑴ 计算数据均值x_mean = np.mean(x)# ⑵ 计算自适应调节因子α# V(x) = 1/n * Σ(xi - x_mean)^2V = np.mean((x - x_mean)**2)alpha = 1 / V# ⑶ 计算权重wi = exp(-(xi - x_mean)^2)weights = np.exp(-alpha * (x - x_mean)**2)# ⑷ 计算位置参数γ = Σ(wi * xi) / Σwigamma = np.sum(weights * x) / np.sum(weights)return gamma
3. 计算形状参数δ
def calculate_delta(x, window_size=5):"""计算形状参数δ参数:x -- 输入数据数组window_size -- 局部计算窗口大小(默认为5)返回:delta -- 形状参数δ"""n = len(x)mu_loc = np.zeros(n)sigma_loc = np.zeros(n)# 使用滑动窗口计算局部均值和标准差for i in range(n):# 确定窗口边界start = max(0, i - window_size // 2)end = min(n, i + window_size // 2 + 1)window_data = x[start:end]# ⑴ 计算局部均值mu_loc[i] = np.mean(window_data)# ⑵ 计算局部标准差if len(window_data) > 1:sigma_loc[i] = np.std(window_data, ddof=1)else:sigma_loc[i] = 0# ⑶ 计算形状参数δ = 局部标准差均值 / 局部均值均值valid_indices = (mu_loc != 0) & (sigma_loc != 0)if np.any(valid_indices):delta = np.mean(sigma_loc[valid_indices]) / np.mean(mu_loc[valid_indices])else:delta = 1.0 # 默认值return delta
4. 计算尺度参数ξ
def calculate_xi(x):"""计算尺度参数ξ参数:x -- 输入数据数组返回:xi -- 尺度参数ξ"""# ⑴ 计算数据均值x_mean = np.mean(x)# ⑵ 计算样本峰度n = len(x)if n < 4:return 1.0 # 样本太小无法计算峰度numerator = np.sum((x - x_mean)**4) / ndenominator = (np.sum((x - x_mean)**2) / n)**2kurtosis = numerator / denominator - 3# ⑶ 计算绝对偏差abs_dev = np.abs(x - x_mean)# ⑷ 计算尺度参数ξ = mean(|xi - x_mean|) * (1 + K(x)/3)/4xi = np.mean(abs_dev) * (1 + kurtosis / 3) / 4return xi
5. 计算伸缩参数λ
def calculate_lambda(x):"""计算伸缩参数λ参数:x -- 输入数据数组返回:lambda_ -- 伸缩参数λ"""# ⑴ 计算数据均值x_mean = np.mean(x)# ⑵ 计算数据标准差sigma = np.std(x, ddof=1)# ⑶ 计算偏态系数if sigma == 0:return 1.0 # 避免除以0skewness = np.sum((x - x_mean)**3) / (n * sigma**3)# ⑷ 计算伸缩参数λ = (1 + |Sk(x)|)/2lambda_ = (1 + np.abs(skewness)) / 2return lambda_
6. Johnson SU 分布概率密度函数
def johnson_su_pdf(x, gamma, delta, xi, lambda_):"""Johnson SU分布概率密度函数参数:x -- 输入数据点gamma -- 位置参数delta -- 形状参数xi -- 尺度参数lambda_ -- 伸缩参数返回:pdf值"""term = (x - gamma) / xiz = gamma + delta * np.arcsinh(term)numerator = delta / (xi * np.sqrt(2 * np.pi))denominator = np.sqrt(1 + term**2)exp_term = np.exp(-0.5 * z**2)return (numerator / denominator) * exp_term
7. 使用坐标下降法优化Johnson SU参数
def negative_log_likelihood(params, x):"""计算负对数似然函数(用于最小化)参数:params -- 参数数组[gamma, delta, xi, lambda_]x -- 观测数据返回:负对数似然值"""gamma, delta, xi, lambda_ = params# 确保参数在合理范围内if delta <= 0 or xi <= 0 or lambda_ <= 0:return np.inf# 计算概率密度pdf_values = johnson_su_pdf(x, gamma, delta, xi, lambda_)# 避免数值问题(取对数前处理0值)pdf_values = np.clip(pdf_values, 1e-10, None)# 计算对数似然log_likelihood = np.sum(np.log(pdf_values))return -log_likelihood
def optimize_johnson_su_params(x, initial_params=None, bounds=None):"""使用坐标下降法优化Johnson SU分布参数参数:x -- 观测数据initial_params -- 初始参数猜测(可选)bounds -- 参数边界(可选)返回:优化后的参数数组[gamma, delta, xi, lambda_]"""# 如果没有提供初始参数，使用我们之前的方法计算if initial_params is None:gamma_init = calculate_gamma(x)delta_init = calculate_delta(x)xi_init = calculate_xi(x)lambda_init = calculate_lambda(x)initial_params = [gamma_init, delta_init, xi_init, lambda_init]# 设置参数边界if bounds is None:# gamma无限制，其他参数必须为正bounds = [(None, None), (1e-6, None), (1e-6, None), (1e-6, None)]# 使用L-BFGS-B算法最小化负对数似然result = minimize(negative_log_likelihood, x0=initial_params, args=(x,),bounds=bounds,method='L-BFGS-B')return result.x
8. 弹性网络参数优化(贝叶斯优化)
def elasticnet_cv(X, y, alpha, l1_ratio):"""弹性网络交叉验证函数参数:X -- 特征矩阵y -- 目标变量alpha -- 正则化强度l1_ratio -- L1正则化比例返回:交叉验证平均得分"""model = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, random_state=42)scores = cross_val_score(model, X, y, cv=5, scoring='neg_mean_squared_error')return np.mean(scores)
def optimize_elasticnet_params(X, y):"""使用贝叶斯优化优化弹性网络参数参数:X -- 特征矩阵y -- 目标变量返回:最佳参数字典"""# 定义参数搜索空间pbounds = {'alpha': (0.0001, 1.0),'l1_ratio': (0.1, 0.9)}# 创建优化器optimizer = BayesianOptimization(f=lambda alpha, l1_ratio: elasticnet_cv(X, y, alpha, l1_ratio),pbounds=pbounds,random_state=42)# 执行优化optimizer.maximize(init_points=5, n_iter=20)return optimizer.max['params']
9. 示例使用
# 生成示例数据
np.random.seed(42)
data = np.random.normal(loc=5, scale=2, size=1000)
# 计算初始参数
gamma = calculate_gamma(data)
delta = calculate_delta(data)
xi = calculate_xi(data)
lambda_ = calculate_lambda(data)
print(f"初始参数: gamma={gamma:.4f}, delta={delta:.4f}, xi={xi:.4f}, lambda={lambda_:.4f}")
# 优化Johnson SU参数
optimized_params = optimize_johnson_su_params(data)
print(f"优化后参数: gamma={optimized_params[0]:.4f}, delta={optimized_params[1]:.4f}, "f"xi={optimized_params[2]:.4f}, lambda={optimized_params[3]:.4f}")
# 弹性网络示例(需要真实数据)
# X = ... # 特征矩阵
# y = ... # 目标变量
# best_params = optimize_elasticnet_params(X, y)
# print("最佳弹性网络参数:", best_params)
代码解释
- 位置参数γ计算:
- 首先计算数据均值
- 然后计算自适应调节因子α作为方差的倒数
- 使用高斯核函数计算每个数据点的权重
- 最后计算加权平均值作为位置参数γ
- 形状参数δ计算:
- 使用滑动窗口计算局部均值和标准差
- 取所有局部标准差与局部均值的比值的平均值作为形状参数δ
- 尺度参数ξ计算:
- 基于数据的峰度和绝对偏差计算
- 峰度衡量数据分布的尾部重量
- 绝对偏差提供数据离散程度的度量
- 伸缩参数λ计算:
- 基于数据的偏态系数计算
- 偏态系数衡量数据分布的不对称性
- Johnson SU分布PDF:
- 实现了Johnson SU分布的概率密度函数公式
- 包含双曲正弦变换和正态分布的组合
- 参数优化:
- 使用负对数似然作为目标函数
- 采用L-BFGS-B算法进行约束优化
- 确保形状、尺度和伸缩参数保持正值
- 弹性网络优化:
- 使用贝叶斯优化寻找最佳正则化参数
- 交叉验证评估模型性能
- 平衡L1和L2正则化的比例
这个实现完整地复现了描述的所有计算步骤，并提供了参数优化方法，可以用于实际数据拟合和分析。