• 基于哈尔小波基的一维密度估计（Python）

    先说点其他的东西。

    关于强非线性、强间断、多物理场强耦合或高度复杂几何形态问题能够得以有效求解的核心难题之一，是如何构建在多尺度情形、非线性作用下具有准确地识别、定位、捕获以及分离各个尺度特征尤其是小尺度局部特征能力的数值工具，这之中包括了大尺度低阶近似解与小尺度高阶微小截断误差的分离解耦。例如，对于湍流等强非线性问题，其核心就在于如何有效地从所关心的大尺度近似解中剥离小尺度特征。对于激波等强间断问题，其核心就在于准确地得到无非物理数值振荡的激波面，且在光滑区域不引入过多的数值粘性。对于多场耦合问题，其难点在于如何以稀疏的数据信息准确表征解和有效地传递局部细节特征。而对于复杂的几何形态，则关键在于有效地刻画出局部几何细节。针对这些问题求解中所体现出的共性需求——即对函数时频局部化特征进行提取与分析的数学能力，小波多分辨分析提供了一种非常有效的数学工具。我们可以通过小波分解系数量级识别局部大梯度特征和局部高频显著信息，且能够通过各种滤波手段得到数值解的稀疏表征，设计出自适应小波数值算法。

    其次，开始基于哈尔小波基的一维密度估计，代码较为简单。

    1. import numpy as np
    2. import matplotlib.pyplot as plt
    3. from scipy.stats import norm, uniform
    4. from scipy.integrate import quad
    5.  
    6.  
    7. import sys, os, time, fileinput
    8.  
    9.  
    10. import matplotlib.pyplot as plt
    11. import matplotlib as mpl
    12.  
    13.  
    14. plt.style.use('default')
    15. # sample data from normal distribution 
    16. N_data = 10000
    17.  
    18.  
    19. # preload sample data
    20. x_data = np.array([])
    21.  
    22.  
    23. # point source samples
    24. Ng = 5
    25. N_scales = uniform.rvs(size = Ng)
    26. scales = 0.005 * uniform.rvs(size = Ng)
    27. scales = 0.005 * np.ones(Ng)
    28. locs = 0.8 * uniform.rvs(size = Ng) + 0.1
    29. for n in range(Ng):
    30.     nm_small = norm(scale = scales[n], loc = locs[n])
    31.     x_data_small = nm_small.rvs(size = int(N_scales[n] * N_data / 20))
    32.     x_data = np.concatenate((x_data, x_data_small))
    33.     
    34. # gaussian samples
    35. nm_large = norm(scale = 0.1, loc = 0.5)
    36. x_data_large = nm_large.rvs(size = N_data)
    37. x_data = np.concatenate((x_data, x_data_large))
    38.  
    39.  
    40. # uniform samples
    41. uni = uniform
    42. x_data_uni = uni.rvs(size = int(N_data / 2))
    43. x_data = np.concatenate((x_data, x_data_uni))
    44.  
    45.  
    46. # plot histogram of distribution
    47. counts, bins, _ = plt.hist(x_data, bins = 100, density = True)
    48. plt.xlim([0,1])
    49.  
    50.  
    51. N = len(x_data)

    1. # load pywt haar wavelets for comparison
    2.  
    3.  
    4. import pywt
    5. wavelet_name = 'haar'
    6.  
    7.  
    8. wavelet = pywt.Wavelet(wavelet_name)
    9. phi, psi, x = wavelet.wavefun(level=9) # level does not affect timing
    10. L = int(x[-1] - x[0])
    11. # rescale to unit length
    12. x = x/L
    13. phi = np.sqrt(L) * phi
    14. psi = np.sqrt(L) * psi
    15.  
    16.  
    17. # define standard haar wavelet and scaling function
    18. def haar_mother_(t): 
    19.     return (np.heaviside(t,0) - np.heaviside(t-1,0) ) * np.interp(t,x,psi)
    20.  
    21.  
    22. def haar_scale_(t):
    23.     return (np.heaviside(t,0) - np.heaviside(t-1,0) ) * np.interp(t,x,phi)
    24.  
    25.  
    26. x_p = np.linspace(-1,1,1000)
    27. plt.plot(x_p, haar_scale_(x_p - 0), c = 'red', alpha = 0.5)
    28. plt.plot(x_p, haar_scale_(2 * x_p - 1), c = 'lightblue', alpha = 0.5)
    29. plt.plot(x_p, haar_scale_(2 * x_p - 0) - haar_scale_(2 * x_p - 1), c = 'purple')
    30.  
    31.  
    32. plt.xlim([-0.25,1.25])

    1. # define dyadic version of wavelet and scaling function
    2. def psi_(x,j,k):
    3.     return 2**(j/2) * haar_mother_(2**j * x  - k)
    4.  
    5.  
    6. def phi_(x,j0,k):
    7.     return 2**(j0/2) * haar_scale_(2**j0 * x  - k)
    8. j1 = 7 # maximum scale (6 -> ~2 min , 7 -> ~9 min)
    9. klist1 = np.arange(0,2**j1 - 1 + 0.5,1) # translations
    10. Nk1 = np.size(klist1)
    11.  
    12.  
    13. scale_info = [j1, klist1]
    14. # plot the density estimate using scaling coefficients
    15.  
    16.  
    17. def angles_for_equal_components_(N):
    18.     """ Generate the angles in spherical coordinates corresponding to a 
    19.     vector whose components are all equal """
    20.     # N = Number of angles required to specify sphere
    21.     arg_for_sqrt = np.arange(N+1, 3-0.5, -1)
    22.     polar_list = np.arccos( 1 /
    23.                            np.sqrt(arg_for_sqrt)
    24.                           )
    25.     azi_list = np.array([np.pi/4])
    26.     return np.concatenate((polar_list, azi_list))
    27.  
    28.  
    29. def initial_scaling_angle_generator_(N, choice):
    30.     """Generates a set of initial scaling angles on the sphere
    31.     N = Dimension of Euclidean Space
    32.     choice = Determine type of scaling angles to generate
    33.     'random': Generate a random sample of angles on the sphere
    34.     'unbiased': All scaling coefficients have the same positive value
    35.     'equiangle': All angles correspond to pi/4"""
    36.     if choice == 'random':
    37.         return np.concatenate((np.pi * np.random.random(size = N - 2), 2*np.pi*np.random.random(size = 1) ))
    38.     elif choice == 'unbiased':
    39.         return angles_for_equal_components_(N-1)
    40.     elif choice == 'equiangle':
    41.         return np.concatenate((np.pi/4 * np.random.random(size = N - 2), np.pi/4*np.ones(1) ))
    42.     
    43. def ct(r, arr):
    44.     """
    45.     coordinate transformation from spherical to cartesian coordinates
    46.     """
    47.     a = np.concatenate((np.array([2*np.pi]), arr))
    48.     si = np.sin(a)
    49.     si[0] = 1
    50.     si = np.cumprod(si)
    51.     co = np.cos(a)
    52.     co = np.roll(co, -1)
    53.     return si*co*r
    54.  
    55.  
    56. def scaling_coefficients_(scaling_angles):
    57.     """
    58.     convert scaling angles in hypersphere to scaling coefficients
    59.     """
    60.     return ct(1,scaling_angles)
    61.  
    62.  
    63. def sqrt_p_(x_data, scaling_angles):
    64.     """
    65.     Calculate the square root of the density estimate as the wavelet expansion
    66.     with the scaling coefficients denoted by the scaling angles
    67.     """
    68.     N = len(x_data)
    69.     j1, klist1 = scale_info
    70.     phi_arr = np.array([phi_(x_data,j1,klist1[nk]) for nk in range(Nk1)])
    71.     scaling_coefficients = scaling_coefficients_(scaling_angles)
    72.     scaling_coefficients_mat = np.outer(scaling_coefficients, np.ones((N,)))
    73.     scale_terms = scaling_coefficients_mat * phi_arr
    74.     return np.sum(scale_terms, axis = 0)
    75.  
    76.  
    77. def safe_log_(x):
    78.     # guard against x = 0
    79.     return np.log(x + 1e-323)
    80.  
    81.  
    82. def unbinned_nll_(x):
    83.     return -np.sum(safe_log_(x))
    84.  
    85.  
    86. def unbinned_nll_sqrt_p_(scaling_angles):
    87.     sqrt_p_arr = sqrt_p_(x_data, scaling_angles)
    88.     p_arr = sqrt_p_arr * sqrt_p_arr
    89.     return unbinned_nll_(p_arr) / N
    90. from iminuit import cost, Minuit
    91.  
    92.  
    93. # unbiased initial scaling angles
    94. initial_scaling_angles = initial_scaling_angle_generator_(Nk1, 'unbiased')
    95. # initial_scaling_angles = initial_scaling_angle_generator_(Nk1, 'random')
    96.  
    97.  
    98. # define iminuit object for minimization
    99. m = Minuit(unbinned_nll_sqrt_p_, initial_scaling_angles)
    100.  
    101.  
    102. # limited to sphere
    103. limits_theta = [(0,np.pi) for n in range(Nk1-2)]
    104. limits_phi = [(0,2*np.pi)]
    105. limits = np.concatenate((limits_theta, limits_phi))
    106.  
    107.  
    108. # set limits
    109. m.limits = limits
    110.  
    111.  
    112. # run fit
    113. m.errordef = Minuit.LIKELIHOOD
    114. m.migrad()

    Migrad
    FCN = -0.3265Nfcn = 10376
    EDM = 8.51e-09 (Goal: 0.0001)time = 544.2 sec
    Valid MinimumBelow EDM threshold (goal x 10)
    SOME parameters at limitBelow call limit
    Hesse okCovariance accurate

    1. # plot the fit
    2. x_plot = np.linspace(0,1,128)
    3. scaling_angles_fit = m.values[:]
    4. sqrt_p_plot = sqrt_p_(x_plot, scaling_angles_fit)
    5. p_plot = sqrt_p_plot * sqrt_p_plot
    6.  
    7.  
    8. plt.plot(x_plot, p_plot, label = 'fit')
    9. counts, bins, _ = plt.hist(x_data, bins = 128, density = True, label = 'data')
    10.  
    11.  
    12. plt.xlim([0,1])
    13. plt.legend()

    1. sqrt_p_plot_j1 = sqrt_p_plot
    2.  
    3.  
    4. j0 = j1 - 4
    5. level = j1 - j0
    6. scaling_coefficients_fit = np.abs(scaling_coefficients_(scaling_angles_fit))
    7. coeffs = pywt.wavedec(scaling_coefficients_fit, wavelet_name, level=level)
    8. scaling_coefficients_j0 = coeffs[0]
    9. wavelet_coefficients_j0 = coeffs[1:]
    10.  
    11.  
    12. klist0 = np.arange(0,2**j0 - 1 + 0.5,1)
    13. Nk0 = np.size(klist0)
    14.  
    15.  
    16. scale_info = [j0, klist0]
    17.  
    18.  
    19. def sqrt_p_coarse_(x_data, scaling_coefficients_j0):
    20.     '''
    21.     Coarse representation of the density estimate
    22.     '''
    23.     N = len(x_data)
    24.     j1, klist1 = scale_info
    25.     phi_arr = np.array([phi_(x_data,j0,klist0[nk]) for nk in range(Nk0)])
    26.     scaling_coefficients_mat = np.outer(scaling_coefficients_j0, np.ones((N,)))
    27.     scale_terms = scaling_coefficients_mat * phi_arr
    28.     return np.sum(scale_terms, axis = 0)
    29.  
    30.  
    31. x_plot = np.linspace(0,1,128)
    32. sqrt_p_plot_j0 = sqrt_p_coarse_(x_plot, scaling_coefficients_j0)
    33. p_plot_j0 = sqrt_p_plot_j0 * sqrt_p_plot_j0
    34.  
    35.  
    36. plt.plot(x_plot, p_plot_j0)
    37. plt.plot(x_plot, p_plot)
    38. counts, bins, _ = plt.hist(x_data, bins = 100, density = True)
    39. plt.xlim([0,1])
    40. plt.legend(["Coarse", "Fine", "Data"])

    1. # Fit summary using DWT
    2.  
    3.  
    4. counts, bins, _ = plt.hist(x_data, bins = 100, density = True)
    5. plt.plot(x_plot, p_plot_j0, label = 'Coarse')
    6. plt.plot(x_plot, p_plot, label = 'Fine')
    7. plt.plot(x_plot, sqrt_p_plot**2 - sqrt_p_plot_j0**2, label = '$p_f - p_c$')
    8. plt.xlim([0,1])
    9. plt.xlabel('$x$')
    10. plt.ylabel('$p(x)$')
    11. plt.legend()

    1. # stationary wavelet transform (test)
    2. # Note: You can't just take the scaling coefficients from the coarse representation; 
    3. ## need to take the inverse (see next cell)
    4. # This leads to a shifted signal (idk why)
    5.  
    6.  
    7. sqrt_p_plot_j1 = sqrt_p_plot
    8.  
    9.  
    10. j0 = j1 - 4
    11. level = j1 - j0
    12. scaling_coefficients_fit = np.abs(scaling_coefficients_(scaling_angles_fit))
    13. coeffs = pywt.swt(scaling_coefficients_fit, wavelet_name, level=level, norm = True, trim_approx = False)
    14. scaling_coefficients_j0 = coeffs[0][0]
    15. wavelet_coefficients_j0 = coeffs[1:]
    16.  
    17.  
    18. klist0 = np.arange(0,2**j1 - 1 + 0.5,1)
    19. Nk0 = np.size(klist1)
    20.  
    21.  
    22. scale_info = [j1, klist1]
    23.  
    24.  
    25. def sqrt_p_coarse_(x_data, scaling_coefficients_j0):
    26.     N = len(x_data)
    27.     j1, klist1 = scale_info
    28.     phi_arr = np.array([phi_(x_data,j1,klist0[nk]) for nk in range(Nk0)])
    29.     scaling_coefficients_mat = np.outer(scaling_coefficients_j0, np.ones((N,)))
    30.     scale_terms = scaling_coefficients_mat * phi_arr
    31.     return np.sum(scale_terms, axis = 0)
    32.  
    33.  
    34. x_plot = np.linspace(0,1,128)
    35. sqrt_p_plot_j0 = sqrt_p_coarse_(x_plot, scaling_coefficients_j0)
    36. p_plot_j0 = sqrt_p_plot_j0 * sqrt_p_plot_j0
    37.  
    38.  
    39. plt.plot(x_plot, p_plot_j0)
    40. plt.plot(x_plot, p_plot)
    41. counts, bins, _ = plt.hist(x_data, bins = 100, density = True)
    42. plt.xlim([0,1])
    43. plt.legend(["Coarse", "Fine", "Data"])

    1. 知乎学术咨询：
    2. https://www.zhihu.com/consult/people/792359672131756032?isMe=1

    工学博士，担任《Mechanical System and Signal Processing》《中国电机工程学报》《控制与决策》等期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。

  • 相关阅读:
    python二次开发CATIA：根据已知数据点创建曲线
    数字藏品NFT交易系统开发 实现元宇宙虚拟事物资产化
    RTX3080安装DEEPLABCUT使用GPU（win11）
    数字图像基础
    Qt 5.12.12 静态编译（MinGW）
    m分别通过matlab和FPGA实现基于高阶循环谱的信号载波调制识别(四阶循环累量)仿真
    swin Transformer
    IC刷卡数据如何获取？5个渠道的IC刷卡数据免费下载~
    Android文件关联
    ubuntu 20.04.4+uWSGI+Nginx安装部署Django+Vue的web前后端全过程记录(1)
  • 原文地址：https://blog.csdn.net/weixin_39402231/article/details/140021919