基于哈尔小波基的一维密度估计(Python)

news2024/11/23 22:03:20

先说点其他的东西。

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

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, uniform
from scipy.integrate import quad


import sys, os, time, fileinput


import matplotlib.pyplot as plt
import matplotlib as mpl


plt.style.use('default')
# sample data from normal distribution 
N_data = 10000


# preload sample data
x_data = np.array([])


# point source samples
Ng = 5
N_scales = uniform.rvs(size = Ng)
scales = 0.005 * uniform.rvs(size = Ng)
scales = 0.005 * np.ones(Ng)
locs = 0.8 * uniform.rvs(size = Ng) + 0.1
for n in range(Ng):
    nm_small = norm(scale = scales[n], loc = locs[n])
    x_data_small = nm_small.rvs(size = int(N_scales[n] * N_data / 20))
    x_data = np.concatenate((x_data, x_data_small))
    
# gaussian samples
nm_large = norm(scale = 0.1, loc = 0.5)
x_data_large = nm_large.rvs(size = N_data)
x_data = np.concatenate((x_data, x_data_large))


# uniform samples
uni = uniform
x_data_uni = uni.rvs(size = int(N_data / 2))
x_data = np.concatenate((x_data, x_data_uni))


# plot histogram of distribution
counts, bins, _ = plt.hist(x_data, bins = 100, density = True)
plt.xlim([0,1])


N = len(x_data)

# load pywt haar wavelets for comparison


import pywt
wavelet_name = 'haar'


wavelet = pywt.Wavelet(wavelet_name)
phi, psi, x = wavelet.wavefun(level=9) # level does not affect timing
L = int(x[-1] - x[0])
# rescale to unit length
x = x/L
phi = np.sqrt(L) * phi
psi = np.sqrt(L) * psi


# define standard haar wavelet and scaling function
def haar_mother_(t): 
    return (np.heaviside(t,0) - np.heaviside(t-1,0) ) * np.interp(t,x,psi)


def haar_scale_(t):
    return (np.heaviside(t,0) - np.heaviside(t-1,0) ) * np.interp(t,x,phi)


x_p = np.linspace(-1,1,1000)
plt.plot(x_p, haar_scale_(x_p - 0), c = 'red', alpha = 0.5)
plt.plot(x_p, haar_scale_(2 * x_p - 1), c = 'lightblue', alpha = 0.5)
plt.plot(x_p, haar_scale_(2 * x_p - 0) - haar_scale_(2 * x_p - 1), c = 'purple')


plt.xlim([-0.25,1.25])

# define dyadic version of wavelet and scaling function
def psi_(x,j,k):
    return 2**(j/2) * haar_mother_(2**j * x  - k)


def phi_(x,j0,k):
    return 2**(j0/2) * haar_scale_(2**j0 * x  - k)
j1 = 7 # maximum scale (6 -> ~2 min , 7 -> ~9 min)
klist1 = np.arange(0,2**j1 - 1 + 0.5,1) # translations
Nk1 = np.size(klist1)


scale_info = [j1, klist1]
# plot the density estimate using scaling coefficients


def angles_for_equal_components_(N):
    """ Generate the angles in spherical coordinates corresponding to a 
    vector whose components are all equal """
    # N = Number of angles required to specify sphere
    arg_for_sqrt = np.arange(N+1, 3-0.5, -1)
    polar_list = np.arccos( 1 /
                           np.sqrt(arg_for_sqrt)
                          )
    azi_list = np.array([np.pi/4])
    return np.concatenate((polar_list, azi_list))


def initial_scaling_angle_generator_(N, choice):
    """Generates a set of initial scaling angles on the sphere
    N = Dimension of Euclidean Space
    choice = Determine type of scaling angles to generate
    'random': Generate a random sample of angles on the sphere
    'unbiased': All scaling coefficients have the same positive value
    'equiangle': All angles correspond to pi/4"""
    if choice == 'random':
        return np.concatenate((np.pi * np.random.random(size = N - 2), 2*np.pi*np.random.random(size = 1) ))
    elif choice == 'unbiased':
        return angles_for_equal_components_(N-1)
    elif choice == 'equiangle':
        return np.concatenate((np.pi/4 * np.random.random(size = N - 2), np.pi/4*np.ones(1) ))
    
def ct(r, arr):
    """
    coordinate transformation from spherical to cartesian coordinates
    """
    a = np.concatenate((np.array([2*np.pi]), arr))
    si = np.sin(a)
    si[0] = 1
    si = np.cumprod(si)
    co = np.cos(a)
    co = np.roll(co, -1)
    return si*co*r


def scaling_coefficients_(scaling_angles):
    """
    convert scaling angles in hypersphere to scaling coefficients
    """
    return ct(1,scaling_angles)


def sqrt_p_(x_data, scaling_angles):
    """
    Calculate the square root of the density estimate as the wavelet expansion
    with the scaling coefficients denoted by the scaling angles
    """
    N = len(x_data)
    j1, klist1 = scale_info
    phi_arr = np.array([phi_(x_data,j1,klist1[nk]) for nk in range(Nk1)])
    scaling_coefficients = scaling_coefficients_(scaling_angles)
    scaling_coefficients_mat = np.outer(scaling_coefficients, np.ones((N,)))
    scale_terms = scaling_coefficients_mat * phi_arr
    return np.sum(scale_terms, axis = 0)


def safe_log_(x):
    # guard against x = 0
    return np.log(x + 1e-323)


def unbinned_nll_(x):
    return -np.sum(safe_log_(x))


def unbinned_nll_sqrt_p_(scaling_angles):
    sqrt_p_arr = sqrt_p_(x_data, scaling_angles)
    p_arr = sqrt_p_arr * sqrt_p_arr
    return unbinned_nll_(p_arr) / N
from iminuit import cost, Minuit


# unbiased initial scaling angles
initial_scaling_angles = initial_scaling_angle_generator_(Nk1, 'unbiased')
# initial_scaling_angles = initial_scaling_angle_generator_(Nk1, 'random')


# define iminuit object for minimization
m = Minuit(unbinned_nll_sqrt_p_, initial_scaling_angles)


# limited to sphere
limits_theta = [(0,np.pi) for n in range(Nk1-2)]
limits_phi = [(0,2*np.pi)]
limits = np.concatenate((limits_theta, limits_phi))


# set limits
m.limits = limits


# run fit
m.errordef = Minuit.LIKELIHOOD
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

# plot the fit
x_plot = np.linspace(0,1,128)
scaling_angles_fit = m.values[:]
sqrt_p_plot = sqrt_p_(x_plot, scaling_angles_fit)
p_plot = sqrt_p_plot * sqrt_p_plot


plt.plot(x_plot, p_plot, label = 'fit')
counts, bins, _ = plt.hist(x_data, bins = 128, density = True, label = 'data')


plt.xlim([0,1])
plt.legend()

sqrt_p_plot_j1 = sqrt_p_plot


j0 = j1 - 4
level = j1 - j0
scaling_coefficients_fit = np.abs(scaling_coefficients_(scaling_angles_fit))
coeffs = pywt.wavedec(scaling_coefficients_fit, wavelet_name, level=level)
scaling_coefficients_j0 = coeffs[0]
wavelet_coefficients_j0 = coeffs[1:]


klist0 = np.arange(0,2**j0 - 1 + 0.5,1)
Nk0 = np.size(klist0)


scale_info = [j0, klist0]


def sqrt_p_coarse_(x_data, scaling_coefficients_j0):
    '''
    Coarse representation of the density estimate
    '''
    N = len(x_data)
    j1, klist1 = scale_info
    phi_arr = np.array([phi_(x_data,j0,klist0[nk]) for nk in range(Nk0)])
    scaling_coefficients_mat = np.outer(scaling_coefficients_j0, np.ones((N,)))
    scale_terms = scaling_coefficients_mat * phi_arr
    return np.sum(scale_terms, axis = 0)


x_plot = np.linspace(0,1,128)
sqrt_p_plot_j0 = sqrt_p_coarse_(x_plot, scaling_coefficients_j0)
p_plot_j0 = sqrt_p_plot_j0 * sqrt_p_plot_j0


plt.plot(x_plot, p_plot_j0)
plt.plot(x_plot, p_plot)
counts, bins, _ = plt.hist(x_data, bins = 100, density = True)
plt.xlim([0,1])
plt.legend(["Coarse", "Fine", "Data"])

# Fit summary using DWT


counts, bins, _ = plt.hist(x_data, bins = 100, density = True)
plt.plot(x_plot, p_plot_j0, label = 'Coarse')
plt.plot(x_plot, p_plot, label = 'Fine')
plt.plot(x_plot, sqrt_p_plot**2 - sqrt_p_plot_j0**2, label = '$p_f - p_c$')
plt.xlim([0,1])
plt.xlabel('$x$')
plt.ylabel('$p(x)$')
plt.legend()

# stationary wavelet transform (test)
# Note: You can't just take the scaling coefficients from the coarse representation; 
## need to take the inverse (see next cell)
# This leads to a shifted signal (idk why)


sqrt_p_plot_j1 = sqrt_p_plot


j0 = j1 - 4
level = j1 - j0
scaling_coefficients_fit = np.abs(scaling_coefficients_(scaling_angles_fit))
coeffs = pywt.swt(scaling_coefficients_fit, wavelet_name, level=level, norm = True, trim_approx = False)
scaling_coefficients_j0 = coeffs[0][0]
wavelet_coefficients_j0 = coeffs[1:]


klist0 = np.arange(0,2**j1 - 1 + 0.5,1)
Nk0 = np.size(klist1)


scale_info = [j1, klist1]


def sqrt_p_coarse_(x_data, scaling_coefficients_j0):
    N = len(x_data)
    j1, klist1 = scale_info
    phi_arr = np.array([phi_(x_data,j1,klist0[nk]) for nk in range(Nk0)])
    scaling_coefficients_mat = np.outer(scaling_coefficients_j0, np.ones((N,)))
    scale_terms = scaling_coefficients_mat * phi_arr
    return np.sum(scale_terms, axis = 0)


x_plot = np.linspace(0,1,128)
sqrt_p_plot_j0 = sqrt_p_coarse_(x_plot, scaling_coefficients_j0)
p_plot_j0 = sqrt_p_plot_j0 * sqrt_p_plot_j0


plt.plot(x_plot, p_plot_j0)
plt.plot(x_plot, p_plot)
counts, bins, _ = plt.hist(x_data, bins = 100, density = True)
plt.xlim([0,1])
plt.legend(["Coarse", "Fine", "Data"])

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

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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/1869888.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

C语言单链表的算法之遍历节点

一:什么是遍历 (1)遍历就是把单链表中的各个节点挨个拿出来,就叫遍历 (2)便利的要点:一是不能遗漏,二是不能重复追求效率 二:如何遍历单链表 (1&#xff0…

element-plus 日期选择添加确定按钮

需求&#xff1a;选择日期后&#xff0c;点击确定按钮关闭面板 思路&#xff1a; 使用shortcuts自定义确定和取消按钮选择日期后使用handleOpen()强制开启面板点击确定后使用handleClose()关闭面板 <template><el-date-pickerref"pickerRef"v-model"…

能量智慧流转:全面升级储能电站的智能网关解决方案

监控系统是电化学储能电站的关键组成部分&#xff0c;储能电站也需要相应的监控系统&#xff0c;通过监控系统对储能设备的状态进行监测&#xff0c;实时感知储能设备的健康状态&#xff0c;控制储能设备的充放电功率和时机等&#xff0c; 一个好的监控系统可以实现储能电站安全…

微软发布Phi-3系列语言模型:手机端的强大AI助手

大模型&#xff08;LLMs&#xff09;在处理复杂任务时展现出的巨大潜力&#xff0c;但却需要庞大的计算资源和存储空间&#xff0c;限制了它们在移动设备等资源受限环境中的应用。微软公司最新发布的Phi-3系列语言模型&#xff0c;以其卓越的性能和小巧的体积&#xff0c;打破了…

[面试题]计算机网络

[面试题]Java【基础】[面试题]Java【虚拟机】[面试题]Java【并发】[面试题]Java【集合】[面试题]MySQL[面试题]Maven[面试题]Spring Boot[面试题]Spring Cloud[面试题]Spring MVC[面试题]Spring[面试题]MyBatis[面试题]Nginx[面试题]缓存[面试题]Redis[面试题]消息队列[面试题]…

AppFlow无代码轻松搭建模型Agent

随着大语言模型发展至今&#xff0c;如何深度开发和使用模型也有了各种各样的答案&#xff0c;在这些答案当中&#xff0c;Agent无疑是一个热点回答。 通过模型也各种插件的组合&#xff0c;可以让你的模型应用具备各种能力&#xff0c;例如&#xff0c;通过天气查询插件机票查…

最新Adobe2024全家桶下载,PS/PR/AE/AI/AU/LR/ID详细安装教程

如大家所熟悉的&#xff0c;Adobe全家桶系列常用的软件有Photoshop&#xff08;PS&#xff09;、Premiere&#xff08;PR&#xff09;、After Effects&#xff08;AE&#xff09;、illustrator&#xff08;AI&#xff09;、Audition&#xff08;AU&#xff09;、Lightroom&…

使用 SwiftUI 为 macOS 创建类似于 App Store Connect 的选择器

文章目录 前言创建选择器组件使用选择器组件总结前言 最近,我一直在为我的应用开发一个全新的界面,它可以让你查看 TestFlight 上所有可用的构建,并允许你将它们添加到测试群组中。 作为这项工作的一部分,我需要创建一个组件,允许用户从特定构建中添加和删除测试群组。我…

商场配电新思维:智能网关驱动的自动化管理系统

在商场配电室监控系统中&#xff0c;主要是以无线网络为载体&#xff0c;目的就是便于对变电站等实时监测与控制。其中&#xff0c;4G配电网关非常关键&#xff0c;可以将配电室系统终端上的信息数据及时上传到服务器&#xff0c;再由服务器下达控制指令到各模块中&#xff0c;…

第六十九:iview 表格汇总怎么拿到传过来的数据,而不是自动累加,需要自定义方法

话不多少&#xff0c;先看官方解释 我这个简单&#xff0c;所以所有说明都在图上了 handleSummary({ columns, data }){console.log(columns, data)let sums {}columns.forEach((item,index)>{const key item.key;console.log("key",item)if(index 0){console.…

003 SpringBoot操作ElasticSearch7.x

文章目录 5.SpringBoot集成ElasticSearch7.x1.添加依赖2.yml配置3.创建文档对象4.继承ElasticsearchRepository5.注入ElasticsearchRestTemplate 6.SpringBoot操作ElasticSearch1.ElasticsearchRestTemplate索引操作2.ElasticsearchRepository文档操作3.ElasticsearchRestTempl…

Mybatis 系列全解(3)——全网免费最细最全,手把手教,学完就可做项目!

Mybatis 系列全解&#xff08;3&#xff09; 1. 多对一处理2. 一对多处理3. 动态SQL3.1 什么是动态SQL3.2 搭建环境3.3 IF3.4 Choose(when,otherwise)3.5 Set3.6 SQL片段3.7 Foreach 4. 缓存4.1 简介4.2 Mybatis 缓存4.3 一级缓存4.4 二级缓存4.5 缓存原理 1. 多对一处理 1&am…

携程任我行有什么用?

眼看一直到十月份都没啥假期了 五一出去玩买了几张携程的卡&#xff0c;想着买景点门票、酒店啥的能有优惠&#xff0c;但最后卡里的钱没用完不说&#xff0c;还有几张压根就没用出去 但是我又不想把卡一直闲置在手里&#xff0c;就怕过期了 最后在收卡云上99.1折出掉了&…

vue3中的图片懒加载指令及全局注册

vue3中的图片懒加载指令及全局注册 最近重新刷了一遍黑马的小兔鲜前端项目&#xff0c;发现有个懒加载的指令之前还没有用过。而且写法相对固定&#xff0c;因此记录一下 首先&#xff0c;懒加载&#xff08;Lazy Loading&#xff09;的作用是延迟加载某些资源或组件&#xf…

Python生成图形验证码

文章目录 安装pillow基本用法生成代码 安装pillow pip install pillow 基本用法 特殊字体文字 如下所示&#xff0c;将下载下来的ttf字体文件放到py文件同一文件夹下 分享一个免费下载字体网站&#xff1a;http://www.webpagepublicity.com/free-fonts.html 我选的字体是Baj…

专题页面设计指南:从构思到实现

如何设计专题页&#xff1f;你有什么想法&#xff1f;专题页的设计主要以发扬产品优势为核心。一个好的专题页可以从不同的角度向用户介绍产品&#xff0c;扩大产品的相关优势&#xff0c;表达产品的优势&#xff0c;让用户在短时间内了解产品。因此&#xff0c;在设计详细信息…

纯css写一个动态圣诞老人

效果预览 在这篇文章中&#xff0c;我们将学习如何使用CSS来创建一个生动的圣诞老人动画。通过CSS的魔力&#xff0c;我们可以让圣诞老人在网页上摇摆&#xff0c;仿佛在向我们招手庆祝圣诞节和新年。 实现思路 实现这个效果的关键在于CSS的keyframes动画规则以及各种CSS属性…

Python 基础:用 json 模块存储和读取数据

目录 一、用 json 存储数据二、用 json 读取数据 遇到看不明白的地方&#xff0c;欢迎在评论中留言呐&#xff0c;一起讨论&#xff0c;一起进步&#xff01; 本文参考&#xff1a;《Python编程&#xff1a;从入门到实践&#xff08;第2版&#xff09;》 用户关闭程序时&#…

step6:改用单例模式

文章目录 文章介绍codemain.cppSerialPort.qmlSerialPortHandler.h 文章介绍 案例MF改为单例模式 参考之前写过的关于单例模式的文章单例模式1、单例模式2 code main.cpp qmlRegisterSingletonType(“com.example.serialport”, 1, 0, “SerialPortHandler”, SerialPortHan…

【Python】易错题 [1]

目录 一、选择&#xff1a; 1.列表的复制​编辑 2.函数 二、填空 一、选择&#xff1a; 1.列表的复制 在Python中&#xff0c;列表是可变的数据类型。当将一个列表赋值给另一个变量时&#xff0c;实际上是将这个变量的引用指向原始列表。&#xff08;指针&#xff09;因此&…