SciPy:半圆上的冯米塞斯分布?

2024-04-07

我试图找出定义包裹在半圆上的冯米塞斯分布的最佳方法(我用它来绘制不同浓度的无方向线)。我目前正在使用 SciPy 的 vonmises.rvs()。本质上,我希望能够输入 pi/2 的平均方向,并将分布截断到两侧不超过 pi/2。

我可以使用截断的正态分布,但我会失去 von-mises 的包裹(假设我想要平均方向为 0)

我在研究映射纤维方向的研究论文中看到了这一点,但我不知道如何实现它(在Python中)。我有点不知道从哪里开始。

如果我的 von Mises 定义为(来自 numpy.von Mises):

np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa))

with:

mu, kappa = 0, 4.0

x = np.linspace(-np.pi, np.pi, num=51)

我该如何改变它以使用围绕半圆的环绕?

有这方面经验的人可以提供一些指导吗?


直接数值逆 CDF 采样很有用,它对于有界域的分布应该非常有效。这是代码示例,构建 PDF 和 CDF 表并使用逆 CDF 方法进行采样。当然可以优化和矢量化

代码、Python 3.8、x64 Windows 10

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

def PDF(x, μ, κ):
    return np.exp(κ*np.cos(x - μ))

N = 201

μ = np.pi/2.0
κ = 4.0

xlo = μ - np.pi/2.0
xhi = μ + np.pi/2.0

# PDF normaliztion

I = integrate.quad(lambda x: PDF(x, μ, κ), xlo, xhi)
print(I)
I = I[0]

x = np.linspace(xlo, xhi, N, dtype=np.float64)
step = (xhi-xlo)/(N-1)

p = PDF(x, μ, κ)/I # PDF table

# making CDF table
c = np.zeros(N, dtype=np.float64)

for k in range(1, N):
    c[k] = integrate.quad(lambda x: PDF(x, μ, κ), xlo, x[k])[0] / I

c[N-1] = 1.0 # so random() in [0...1) range would work right

#%%
# sampling from tabular CDF via insverse CDF method

def InvCDFsample(c, x, gen):
    r = gen.random()
    i = np.searchsorted(c, r, side='right')
    q = (r - c[i-1]) / (c[i] - c[i-1])
    return (1.0 - q) * x[i-1] + q * x[i]

# sampling test
RNG = np.random.default_rng()

s = np.empty(20000)

for k in range(0, len(s)):
    s[k] = InvCDFsample(c, x, RNG)

# plotting PDF, CDF and sampling density
plt.plot(x, p, 'b^') # PDF
plt.plot(x, c, 'r.') # CDF
n, bins, patches = plt.hist(s, x, density = True, color ='green', alpha = 0.7)
plt.show()

以及带有 PDF、CDF 和采样直方图的图表

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

SciPy:半圆上的冯米塞斯分布? 的相关文章

  • 如何创建一个在给定范围内随机打乱数字的 int 数组[重复]

    这个问题在这里已经有答案了 基本上 假设我有一个可以容纳 10 个数字的 int 数组 这意味着我可以在每个索引中存储 0 9 每个数字只能存储一次 如果我运行下面的代码 int num new int 10 for int i 0 i l
  • scipy-optimize-minimize 不执行优化 - CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL

    我试图最小化定义如下的函数 utility decision decision risk cost 其中变量采用以下形式 决策 二进制数组 风险 浮点数数组 成本 常数 我知道解决方案将采取以下形式 决定 1如果 风险 gt 阈值 决定 0
  • 将区间映射到更小的区间的算法

    我尝试搜索 但由于问题的性质 我无法找到满意的内容 我的问题如下 我试图将 0 到 2000 范围内的数字 尽管理想情况下上限是可调的 映射到 10 到 100 范围内的更小的区间 上限将映射 2000 gt 100 和下限也是如此 除此之
  • 可重复的随机数系列

    如何在 PHP 中获得一系列可重复的伪随机数 在旧版本的 PHP 中 我只需在RNG http en wikipedia org wiki Random number generation 但它不再起作用了 因为 PHP 改变了 rand
  • 随机数生成器每次仅返回一个数字

    Python 是否有一个随机数生成器 每次只返回一个随机整数next 函数被调用 数字不应该重复并且生成器应返回区间内的随机整数 1 1 000 000 这是独一无二的 我需要生成超过一百万个不同的数字 这听起来好像非常消耗内存 以防所有数
  • perl生成字符串来匹配正则表达式

    我尝试找到一种方法来生成与正则表达式匹配的字符串 例如以下正则表达式 A Z 6 6 A Z2 9 A NP Z0 9 A Z0 9 3 3 0 1 我尝试过 Cpan 上的一些 perl 模块不起作用 gt 字符串 随机 gt 正则表达式
  • 如何使用 Misc.imread 将图像分割为红色、绿色和蓝色通道

    我正在尝试将图像切片为 RGB 但在绘制这些图像时遇到问题 我使用此函数从某个文件夹获取所有图像 def get images path image type image list for filename in glob glob pat
  • 无需时间即可生成随机字符串?

    我知道如何使用 Runes 和播种 rand Init 在 go 中生成随机字符串time UnixNano 我的问题是 是否可以 使用 stdlib 在不使用当前时间戳 安全 的情况下播种 rand 此外 我问 因为仅仅依靠时间来为敏感操
  • 导入 scipy.stats 时,出现“ImportError: DLL load failed: 找不到指定的过程”

    我无法导入 scipy stats 并收到以下错误 但不知何故 import scipy as sp 仍然可以正常工作 其他库如numpy pandas都可以毫无问题地导入 我尝试在 Anaconda 中重新安装 scipy 1 2 1 降
  • `numpy.diff` 和 `scipy.fftpack.diff` 在微分时给出不同的结果

    我正在尝试计算一些数据的导数 并且正在尝试比较有限差分的输出和谱方法的输出 但结果却截然不同 我无法弄清楚到底为什么 考虑下面的示例代码 import numpy as np from scipy import fftpack as sp
  • 在哪里可以获得几乎所有英语单词的列表? [关闭]

    Closed 这个问题不符合堆栈溢出指南 help closed questions 目前不接受答案 我想生成一些随机文本 我尝试写一个基本的Java程序 int nowords r nextInt 2000 int i j for i 0
  • Python 或 C 语言中的 Matlab / Octave bwdist()

    有谁知道 Matlab Octave bwdist 函数的 Python 替代品 此函数返回给定矩阵的每个单元格到最近的非零单元格的欧几里得距离 我看到了一个 Octave C 实现 一个纯 Matlab 实现 我想知道是否有人必须用 AN
  • 如何在 python 中将最佳概率分布模型拟合到我的数据?

    我有大约 20 000 行这样的数据 Id value 1 30 2 3 3 22 n 27 我对我的数据进行了统计 平均值33 85 中位数30 99 最小值2 8 最大值206 95 置信区间0 21 所以大多数值在33左右 并且有一些
  • 如何在Python中计算输出的均值、众数、方差、标准差等?

    我有一个基于概率的简单游戏 每天我们抛一枚硬币 如果正面朝上 我们就赢了 我们会得到 20 美元 如果我们抛硬币 反面朝上 那么我们会在月底损失 19 美元 28 天 我们可以看到我们失去或赚了多少 def coin tossing gam
  • 线性同余生成器 - 如何选择种子和统计检验

    我需要做一个线性同余生成器 它将成功通过所选的统计测试 我的问题是 如何正确选择发电机的数字以及 我应该选择哪些统计检验 我想 均匀性的卡方频率测试 每代收集10 000个号码的方法 将 0 1 细分为10个相等的细分 柯尔莫哥洛夫 斯米尔
  • 沿轴 0 重复 scipy csr 稀疏矩阵

    我想重复 scipy csr 稀疏矩阵的行 但是当我尝试调用 numpy 的重复方法时 它只是将稀疏矩阵视为对象 并且只会将其作为 ndarray 中的对象重复 我浏览了文档 但找不到任何实用程序来重复 scipy csr 稀疏矩阵的行 我
  • 如何指定聚类的距离函数?

    我想对给定距离的点进行聚类 奇怪的是 似乎 scipy 和 sklearn 聚类方法都不允许指定距离函数 例如 在sklearn cluster AgglomerativeClustering 我唯一可以做的就是输入一个亲和力矩阵 这将非常
  • wavfile.read python 文件意外结束

    我正在尝试通过以下代码读取 wav 音频文件 from scipy io import wavfile file PC1 20090513 050000 0010 wav rate audio wavfile read file 但它显示以
  • 使用时间、日期、时间增量

    我有一个问题 我的工作时间和时差很多 到目前为止 我已经使用许多 if 语句解决了这个问题 但这些语句很容易出错 在寻找更好的解决方案并且无需重新发明轮子的过程中 我遇到了时间 日期和时间增量 但这些对我来说似乎太不灵活了 所以我正在寻找如
  • 计算互相关函数?

    In R 我在用ccf or acf计算成对互相关函数 以便我可以找出哪个移位给我带来最大值 从它的外观来看 R给我一个标准化的值序列 Python 的 scipy 中是否有类似的东西 或者我应该使用fft模块 目前 我正在这样做 xcor

随机推荐

  • 将字节数组转换为 CGImage

    我在无符号字符数组中有图像的灰度值 我想将其转换为CGImage这样我就可以使用它在 iOS 中显示它UIImage unsigned char 数组的每个值都是灰度图像的像素值 我正在使用以下代码 但图像未显示 void viewDidL
  • 与不使用 if 的测试相比,if 语句的效率如何? (C++)

    我需要一个程序来获取两个数字中较小的一个 我想知道是否使用标准 如果 x 小于 y int a b low if a lt b low a else low b 比这个更有效率或更低 int a b low low b a b a b gt
  • Excel:找出单元格是否包含/包含一系列单元格中的值

    我有一个职位列表 A 和短语列表 B 对于 A 中的每个标题 我想检查它是否包含 B 中的短语 任何短语 我不在乎是哪个 1 Example Column A Example Column B 2 Head of Marketing Sen
  • Kafka 无法在 Azure Web 应用服务上运行

    我们正在尝试在 Azure Web 应用服务上部署 Kafka 我们不断收到以下错误 容器 NN 未响应端口 2181 上的 HTTP ping 站点启动失败 请参阅容器日志进行调试 WEBSITES PORT 设置为 2181 并且在 D
  • 密码自加密可行吗?

    我一直在阅读有关密码存储的内容 基本上发现了两种常用的技术 使用单个密钥对存储的所有密码进行加密 Using hashes 使用带盐的哈希值 存储 自行加密 的密码是否存在缺陷 即加密一个txt 其中写着password1用密码passwo
  • AngularJS 如何发送多部分/混合

    我正在开发一个项目 我必须在 AngularJS 中上传少量 JSON 和文件 我已经使用 Danial Farid 的 Angular file upload 编写了代码 并且它正在工作 除了它总是发送 multipart form da
  • 运行 flutter doctor 时 Android 许可证状态未知

    我无法运行 flutter 应用程序 因为 Android 许可证未知的 cmd 部分中不断弹出错误 我什至尝试过更新 android studio 但没有帮助 还出现了一个错误 但我交叉检查了它要求删除的文件已经被删除 并且 androi
  • ASP.NET MVC - 使用相同的表单来创建和编辑

    创建用于创建新模型和编辑现有模型的表单的最佳实践方法是什么 人们可以为我指明方向吗 书呆子晚餐 http nerddinner codeplex com will really指明道路 创建 aspx
  • 如何从 webServiceTemplate 获取肥皂响应

    我需要使用 webServiceTemplate 获得肥皂响应 目前 在我现有的架构中 它使用函数 public boolean sendSourceAndReceiveToResult String uri Source requestP
  • 抗锯齿模式差异?

    这两种抗锯齿模式有区别吗 e Graphics SmoothingMode Drawing2D SmoothingMode AntiAlias e Graphics SmoothingMode Drawing2D SmoothingMode
  • Jekyll 中 localhost 和 github 页面的 Baseurl 行为不同

    我正在开发一个静态网站Jekyll 部署于github pages 我在使用配置文件中的 baseurl 时遇到问题 这是我的摘录 config yml baseurl blog url http remidoolaeghe github
  • Python:getopt、批处理文件和带空格的路径

    我正在使用 getopt 来解析选项和参数 我编写了一个批处理文件来调用 python 脚本 这样我就不必一遍又一遍地输入相同的命令 当我打印出参数列表时 路径被空格分割 并且每个路径都被单独解析 该路径用双引号引起来 但我不确定问题是什么
  • 无法配置 HTTPS 端点。未指定服务器证书,找不到默认的开发者证书

    我正在开发一个已配置 HTTPS 的结构应用程序 尽管我有有效的安装证书 但它抛出异常 这些说明来自这个博客 https www waynethompson com au blog dotnet dev certs https 为我工作 d
  • PHP:如何在字符串中的随机位置添加随机字符

    如何在字符串中的随机位置添加单个随机字符 0 9 或 a z 或 或 我可以通过以下方式获得随机位置 random position rand 0 5 现在我怎样才能得到一个随机数 0到9 OR随机字符 a 到 z OR OR 最后 我如何
  • 如何在 Xamarin for Android 中压缩文件?

    我有一个函数可以创建一个 zip 文件和传递的文件字符串数组 该函数确实成功创建了 zip 文件及其内部的 zip 条目文件 但这些 zip 条目文件是空的 我尝试了几种不同的方法 下面的函数代码是我最接近的工作代码 public stat
  • miniconda 和 miniforge 之间有什么区别?

    The 小型锻造厂 https github com conda forge miniforgeinstaller 是一个相对较新的 社区主导的 最小的 conda 安装程序 正如其自述文件中所述 可以直接与 Miniconda 进行比较
  • 如何在 Q/KDB 中生成格式化的日期字符串?

    如何从 Q 日期类型生成 ISO 日期字符串 yyyy MM dd 我考虑过连接各个部分 但我什至无法获取日期 月份 例如d 2015 12 01 d month prints 2015 12 即不仅仅是月份 如果您计划大规模执行此操作 即
  • pip、代理身份验证和“不支持代理方案”

    尝试在新的 python 安装上安装 pip 我遇到了代理错误 看起来像一个错误get pip or urllib3 问题是我是否必须经历设置的痛苦CNTLM 如此处所述 https stackoverflow com questions
  • 如何从控制台输入填充整数数组

    假设我知道用户将输入多少个数字 我有一个 int 数组 我想用用户输入的按特定字符 例如空格 分隔的整数来填充它 我设法用这种方法解决了它 int numbers new int 5 string input Console ReadLin
  • SciPy:半圆上的冯米塞斯分布?

    我试图找出定义包裹在半圆上的冯米塞斯分布的最佳方法 我用它来绘制不同浓度的无方向线 我目前正在使用 SciPy 的 vonmises rvs 本质上 我希望能够输入 pi 2 的平均方向 并将分布截断到两侧不超过 pi 2 我可以使用截断的