Modified Bessel Function of the First Kind

Abstract

最近接触到 von Mises–Fisher distribution, 其概率密度如下:
f p ( x ; μ , κ ) = κ p 2 − 1 ( 2 π ) p 2 I p 2 − 1 ( κ ) e κ μ ⊺ x \begin{aligned} f_{p}(\bm{x}; \bm{\mu}, \kappa) = \frac{\kappa^{\frac{p}{2}-1}} {(2\pi)^{\frac{p}{2}} I_{\frac{p}{2}-1}(\kappa)} e^{\kappa\bm{\mu^\intercal}\bm{x}} \end{aligned} fp(x;μ,κ)=(2π)2pI2p1(κ)κ2p1eκμx 其中 I α ( x ) I_{\alpha}(x) Iα(x) α \alpha α 阶第一类贝塞尔函数.查阅 Wikipedia, 其公式为:


前者是定义式, 后者是积分式. 其计算比较麻烦, 更不用说求导了, 好在这类带有 e ∗ ∗ ∗ e^{***} e∗∗∗ 的积分式能推出递推式, 主要有:

1. 推导 d I α d x = I α − 1 ( x ) − α x I α ( x ) \frac{dI_{\alpha}}{dx} = I_{\alpha-1}(x) - \frac{\alpha}{x} I_{\alpha}(x) dxdIα=Iα1(x)xαIα(x):

I α ( x ) = ∑ m = 0 ∞ ( x 2 ) 2 m + α m ! Γ ( m + α + 1 ) d I α d x = ∑ m = 0 ∞ ( 2 m + α ) ( x 2 ) 2 m + α − 1 2 m ! Γ ( m + α + 1 ) = ∑ m = 0 ∞ ( 2 m + 2 α − α ) ( x 2 ) 2 m + α − 1 2 m ! Γ ( m + α + 1 ) = ∑ m = 0 ∞ ( m + α ) ( x 2 ) 2 m + [ α − 1 ] m ! ( m + α ) Γ ( m + [ α − 1 ] + 1 ) − α x ∑ m = 0 ∞ ( x 2 ) ( x 2 ) 2 m + α − 1 m ! Γ ( m + α + 1 ) = I α − 1 ( x ) − α x I α ( x ) \begin{aligned} I_{\alpha}(x) &= \sum_{m=0}^{\infin} \frac{(\frac{x}{2})^{2m+\alpha}}{m!\Gamma(m+\alpha+1)} \\ \frac{dI_{\alpha}}{dx} &= \sum_{m=0}^{\infin} \frac{(2m+\alpha)(\frac{x}{2})^{2m+\alpha-1}}{2m!\Gamma(m+\alpha+1)} \\ &= \sum_{m=0}^{\infin} \frac{(2m+2\alpha-\alpha)(\frac{x}{2})^{2m+\alpha-1}}{2m!\Gamma(m+\alpha+1)} \\ &= \sum_{m=0}^{\infin} \frac{(m+\alpha)(\frac{x}{2})^{2m+[\alpha-1]}}{m!(m+\alpha)\Gamma(m+[\alpha-1]+1)} - \frac{\alpha}{x}\sum_{m=0}^{\infin} \frac{(\frac{x}{2})(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)} \\ &= I_{\alpha-1}(x) - \frac{\alpha}{x} I_{\alpha}(x) \end{aligned} Iα(x)dxdIα=m=0m!Γ(m+α+1)(2x)2m+α=m=02m!Γ(m+α+1)(2m+α)(2x)2m+α1=m=02m!Γ(m+α+1)(2m+2αα)(2x)2m+α1=m=0m!(m+α)Γ(m+[α1]+1)(m+α)(2x)2m+[α1]xαm=0m!Γ(m+α+1)(2x)(2x)2m+α1=Iα1(x)xαIα(x)

2. 推导 d I α d x = I α + 1 ( x ) + α x I α ( x ) \frac{dI_{\alpha}}{dx} = I_{\alpha+1}(x) + \frac{\alpha}{x} I_{\alpha}(x) dxdIα=Iα+1(x)+xαIα(x):

d I α d x = ∑ m = 0 ∞ ( 2 m + α ) ( x 2 ) 2 m + α − 1 2 m ! Γ ( m + α + 1 ) = ∑ m = 0 ∞ m ( x 2 ) 2 m + α − 1 m ! Γ ( m + α + 1 ) + ∑ m = 0 ∞ α ( x 2 ) 2 m + α − 1 2 m ! Γ ( m + α + 1 ) \begin{aligned} \frac{dI_{\alpha}}{dx} &= \sum_{m=0}^{\infin} \frac{(2m+\alpha)(\frac{x}{2})^{2m+\alpha-1}}{2m!\Gamma(m+\alpha+1)} \\ &= \sum_{m=0}^{\infin} \frac{m(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)} + \sum_{m=0}^{\infin} \frac{\alpha(\frac{x}{2})^{2m+\alpha-1}}{2m!\Gamma(m+\alpha+1)} \\ \end{aligned} dxdIα=m=02m!Γ(m+α+1)(2m+α)(2x)2m+α1=m=0m!Γ(m+α+1)m(2x)2m+α1+m=02m!Γ(m+α+1)α(2x)2m+α1 后一项就是 α x I α ( x ) \frac{\alpha}{x} I_{\alpha}(x) xαIα(x), 但前面一项看着不好搞, 式中 x 2 \frac{x}{2} 2x 的次数是 ( 2 m + α − 1 ) (2m+\alpha-1) (2m+α1), 却能推到 I α + 1 ( x ) I_{\alpha+1}(x) Iα+1(x), 起初还从 Γ ( m + α + 1 ) = Γ ( m + [ α + 1 ] + 1 ) m + α + 1 \Gamma(m+\alpha+1) = \frac{\Gamma(m+[\alpha+1]+1)}{m+\alpha+1} Γ(m+α+1)=m+α+1Γ(m+[α+1]+1) 搞, 但发现结果是: ∑ m = 0 ∞ m ( x 2 ) 2 m + α − 1 m ! Γ ( m + α + 1 ) = ∑ m = 0 ∞ m ( m + α + 1 ) ( x 2 ) 2 m + α − 1 m ! Γ ( m + [ α + 1 ] + 1 ) \sum_{m=0}^{\infin} \frac{m(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)} = \sum_{m=0}^{\infin} \frac{m(m+\alpha+1)(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+[\alpha+1]+1)} m=0m!Γ(m+α+1)m(2x)2m+α1=m=0m!Γ(m+[α+1]+1)m(m+α+1)(2x)2m+α1 次数都不对, 咋搞? 发现若约去 m m m, 则 ( m − 1 ) ! ∣ m = 0 = ( − 1 ) ! (m-1)!|_{m=0} = (-1)! (m1)!m=0=(1)!, 那就单独拿出来看看: m ( x 2 ) 2 m + α − 1 m ! Γ ( m + α + 1 ) ∣ m = 0 = 0 \frac{m(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)}|_{m=0} = 0 m!Γ(m+α+1)m(2x)2m+α1m=0=0, 则:
d I α d x = ∑ m = 0 ∞ m ( x 2 ) 2 m + α − 1 m ! Γ ( m + α + 1 ) + α x I α ( x ) = ∑ m = 1 ∞ ( x 2 ) 2 m + α − 1 ( m − 1 ) ! Γ ( m + α + 1 ) + α x I α ( x ) t = m − 1 ⇒ = ∑ t = 0 ∞ ( x 2 ) 2 ( t + 1 ) + α − 1 t ! Γ ( ( t + 1 ) + α + 1 ) + α x I α ( x ) = ∑ t = 0 ∞ ( x 2 ) 2 t + [ α + 1 ] t ! Γ ( t + [ α + 1 ] + 1 ) + α x I α ( x ) = I α + 1 ( x ) + α x I α ( x ) \begin{aligned} \frac{dI_{\alpha}}{dx} &= \sum_{m=0}^{\infin} \frac{m(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)} + \frac{\alpha}{x} I_{\alpha}(x) \\ &= \sum_{m=1}^{\infin} \frac{(\frac{x}{2})^{2m+\alpha-1}}{(m-1)!\Gamma(m+\alpha+1)} + \frac{\alpha}{x} I_{\alpha}(x) \\ t = m-1 \Rightarrow &= \sum_{t=0}^{\infin} \frac{(\frac{x}{2})^{2(t+1)+\alpha-1}}{t!\Gamma((t+1)+\alpha+1)} + \frac{\alpha}{x} I_{\alpha}(x) \\ &= \sum_{t=0}^{\infin} \frac{(\frac{x}{2})^{2t+[\alpha+1]}}{t!\Gamma(t+[\alpha+1]+1)} + \frac{\alpha}{x} I_{\alpha}(x) \\ &= I_{\alpha+1}(x) + \frac{\alpha}{x} I_{\alpha}(x) \end{aligned} dxdIαt=m1=m=0m!Γ(m+α+1)m(2x)2m+α1+xαIα(x)=m=1(m1)!Γ(m+α+1)(2x)2m+α1+xαIα(x)=t=0t!Γ((t+1)+α+1)(2x)2(t+1)+α1+xαIα(x)=t=0t!Γ(t+[α+1]+1)(2x)2t+[α+1]+xαIα(x)=Iα+1(x)+xαIα(x)

3. 其他两个递推式

推了这两个, 那么两式相加得: I α − 1 ( x ) + I α + 1 ( x ) = 2 d I α d x I_{\alpha-1}(x) + I_{\alpha+1}(x) = 2 \frac{dI_{\alpha}}{dx} Iα1(x)+Iα+1(x)=2dxdIα 两式相减得: I α − 1 ( x ) − I α + 1 ( x ) = 2 I α I_{\alpha-1}(x) - I_{\alpha+1}(x) = 2 I_{\alpha} Iα1(x)Iα+1(x)=2Iα 就很自然了.

4. 积分式的递推式证明

见第一类修正贝塞尔函数的递推公式怎么证明?, 这是阶数 α \alpha α 为整数时的情况.

5. Exponentially Scaled Modified Bessel Function of the First Kind.

Defined as::
	ive(v, z) = iv(v, z) * exp(-abs(z.real))

所以当假定 z . r e a l ≥ 0 z.real \ge 0 z.real0 时, 求导数就按照莱布尼茨法则, 就可以了:

ive(v, z)' = ive(ctx.v - 1, z) - ive(ctx.v, z) * (ctx.v + z) / z

6. 用 scipy.special 包进行验证

import numpy as np
from scipy import special

v = 5.5
z = 3
iv = special.iv(v, z)
ive = special.ive(v, z)  # iv * np.exp(-abs(Re(z)))
print(iv)  # 0.04532335799965591
print(ive)  # 0.0022565171233900425

print(iv * np.exp(-z))  # 验证 ive = iv * exp(-z), 0.002256517123390042

# %% 验证导数
ivp = special.ivp(v, z)
ivp_1 = (special.iv(v - 1, z) + special.iv(v + 1, z)) / 2
ivp_2 = special.iv(v - 1, z) - v / z * iv
ivp_3 = special.iv(v + 1, z) + v / z * iv

print(ivp_1)  # 0.09310529508597687
print(ivp_2)  # 0.09310529508597686
print(ivp_3)  # 0.09310529508597688

# %% 验证ive的导数
ivep = (ivp - iv) * np.exp(-z)
ivep_1 = special.ive(v - 1, z) - (z + v) / z * ive
print(ivep)  # 0.0023789225684656356
print(ivep_1)  # 0.002378922568465644
补充材料

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

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

2.5K Star,打造个性化博客平台

2.5K Star,打造个性化博客平台 Hi,骚年,我是大 G,公众号「GitHub 指北」会推荐 GitHub 上有趣有用的项目,一分钟 get 一个优秀的开源项目,挖掘开源的价值,欢迎关注。 导语 在当今的信息时代&a…

Interceptor拦截器+JWT令牌实现登陆验证

一、背景 与过滤器的作用类似,不过拦截器是spring中的组件,只能拦截进入spring的请求;过滤器则可以拦截所有从前端页面发送来的请求。 *拦截器和过滤器选一就可以实现登陆验证,过滤器的实现在以下这篇博客中,有需要可…

设置Matlab2022a断点查看参数变化

Matlab2022a设置断点,查看参数变化 本文使用的是下载好的matlab2022a软件,下载加安装matlab预计1小时(百度网盘加速)。需要的安装包的评论。 安装好的matlab界面如下: 接下来,编辑一个.m文件,…

【docker基础学习之】镜像构建

下面是在工作过遇到的一些实际例子,谨以此作为笔记参考 目录 1.背景2. 寻找方案3. 如何解决4.解决步骤4.1 DockerFile4.2 现在要做的 5. 镜像相关命令 1.背景 部署(迁移)项目时发现,项目的excel导出功能报错,错误如下…

灵根孕育源流出,心性修持大道生

解法&#xff1a; 手动本地跑了一下1e9&#xff0c;显然超时。 然后预处理发现开不了这么大的数组。 肯定有规律&#xff0c;打表看看 代码如下 #include<iostream> #include<vector> #include<algorithm> #include<cmath> using namespace std; #…

3.7作业

一 1&#xff09;应用层 负责处理不同应用程序之间的通信&#xff0c;需要满足提供的协议&#xff0c;确保数据发送方和接收方的正确 应用层提供的协议&#xff1a; &#xff08;2&#xff09;表示层 负责网络中通信的数据的编码和格式&#xff0c;确保通信过程中…

喜讯|智安网络与OceanBase完成产品兼容性互认证

近日&#xff0c;深圳市智安网络有限公司旗下产品智安云综合防御平台V3.0与北京奥星贝斯科技有限公司的OceanBase数据库正式完成兼容性互认证。经双方联合测试&#xff0c;结果表明&#xff1a;OceanBase数据库软件V4与智安云综合防御平台V3.0完全兼容&#xff0c;在功能、性能…

用户角色的重要性:确保财务数据安全的最佳方式

在企业的财务管理业务中&#xff0c;一个人几乎不可能完成所有的财务记账任务&#xff0c;例如设定预算、发票审批等等&#xff0c;至少不能有效地执行。最为明智的方式&#xff0c;是将这些任务分派给特定的人员&#xff0c;比如部门经理、财务经理或者销售、市场人员等等。 但…

Sora: 大型视觉模型背景、技术、局限性和机遇的综述

论文链接&#xff1a;https://arxiv.org/pdf/2402.17177.pdf 背景 在分析 Sora 之前&#xff0c;研究者首先盘点了视觉内容生成技术的沿袭。 在深度学习革命之前&#xff0c;传统的图像生成技术依赖于基于手工创建特征的纹理合成和纹理映射等方法。这些方法在生成复杂而生动…

Xilinx 7系列 FPGA硬件知识系列(八)——Xilinx FPGA的复位

目录 概要 Xilinx复位准则 全局复位主要由以下三种方式实现 高时钟频率下的复位时序全局复位对时序的要求真的很关键吗&#xff1f; 独热码状态机的复位 FPGA配置 概要 Xilinx白皮书WP272《Get Smart About Reset: Think Local, Not Global》详细讲述了FPGA的全…

IDEA中新增文件,弹出框提示是否添加到Git点错了,怎么重新设置?

打开一个配置了Git的项目&#xff0c;新增一个文件&#xff0c;会弹出下面这个框。提示是否将新增的文件交给Git管理。 一般来说&#xff0c;会选择ADD&#xff0c;并勾选Dont ask agin&#xff0c;添加并不再询问。如果不小心点错了&#xff0c;可在IDEA中重新设置&#xff08…

蓝桥杯嵌入式模板构建——RCT时钟

在CubeMX里的RTC模块启用RTC时钟和日历功能 输入到RTC的时钟要配置成1HZ,这样的话RTC每经过1s走时一次 由于RTC时钟默认配置为32Khz 所以我们需要将异步分频值与同步分频值的乘积调整为32K分频即可一秒走时一次 频率&#xff1a;32000hz / 32000hz 1hz 必须是31和999&#…

Processing基本形状内容和实例

一、Processing的基本形状内容和实例 1.Processing有一组专门绘制基本图形得图案。像线条这样的基本图形可以被连接起来创建更为复杂得形状&#xff0c;例如一片叶子或者一张脸。 2.为了绘制一条直线&#xff0c;我们需要四个参数&#xff0c;两个用于确定初始位置&#xff0c;…

一文读懂HDMI的演变-从HDMI1.0到HDMI2.1(建议收藏)

HDMI&#xff0c;全称为&#xff08;High Definition Multimedia Interface&#xff09;高清多媒体接口&#xff0c;主要用于传输高清音视频信号。 HDMI System HDMI系统包括HDMI的source和HDMI的sink, 其中source 是源端&#xff0c;即信号的来源&#xff1b;Sink的接收端&a…

Android车载开发之AAOS快速入门

一、概述 在正式介绍Android Automotive OS之前,我们先弄清两个概念:Android Auto和Android Automotive OS。 Android Auto Android Auto 不是操作系统,而是一个应用或一个服务。当 Android 手机通过无线或有线方式连接到汽车时,Android 系统会将使用 Android Auto 服务…

python并发编程:阻塞IO

阻塞IO&#xff08;blocking IO&#xff09; 在Linux中&#xff0c;默认情况下所有的socket都是blocking&#xff0c;一个典型的读操作流程大概是这样&#xff1a; 当用户进程调用了recvfrom这个系统调用&#xff0c;kernel就开始了IO的第一个阶段&#xff1a;准备数据。对于…

力扣513 找树左下角的值 Java版本

文章目录 题目描述解题思路代码 题目描述 给定一个二叉树的 根节点 root&#xff0c;请找出该二叉树的 最底层 最左边 节点的值。 假设二叉树中至少有一个节点。 示例 1: 输入: root [2,1,3] 输出: 1 示例 2: 输入: [1,2,3,4,null,5,6,null,null,7] 输出: 7 提示: 二…

Openwrt(IstoreOS)安装iventoy

背景 目前家里有两台不用的旧主机&#xff0c;平时没事在家里折腾这两台机器。经常换装各种系统。最早是将镜像刷入u盘作为启动盘&#xff0c;这样需要重复装系统就特别麻烦。后来用了ventoy以后一个U盘可以放多个系统镜像&#xff0c;还能做口袋系统&#xff08;SystemToGo&a…

MedSAM 项目排坑记录

MedSAM 项目排坑记录 任务排坑过程配置python环境测试构建docker模型训练数据预处理 单GPU训练最后推理 任务 做一个课程大作业&#xff0c;需要进行CVPR2024年医疗影像分割赛题的打榜&#xff08;CVPR 2024: SEGMENT ANYTHING IN MEDICAL IMAGES ON LAPTOP&#xff09;。看到…

Flink实时数仓同步:切片表实战详解

一、背景 在大数据领域&#xff0c;初始阶段业务数据通常被存储于关系型数据库&#xff0c;如MySQL。然而&#xff0c;为满足日常分析和报表等需求&#xff0c;大数据平台采用多种同步方式&#xff0c;以适应这些业务数据的不同存储需求。 一项常见需求是&#xff0c;业务使用…