拟合 Von Mises 分布的圆形直方图

Question

在过去的几天里，我一直在尝试使用 python 绘制圆形数据，方法是构建一个范围从 0 到 2pi 的圆形直方图并拟合 Von Mises 分布。我真正想要实现的是：

1. 具有拟合 Von-Mises 分布的方向数据。该图是使用 Matplotlib、Scipy 和 Numpy 构建的，可以在以下位置找到：http ://jpktd.blogspot.com/2012/11/polar-histogram.html

2. 该图是使用 R 生成的，但给出了我想要绘制的内容的想法。可以在这里找到：https : //www.zeileis.org/news/circtree/

到目前为止我做了什么：

from scipy.special import i0  
import numpy as np
import matplotlib.pyploy as plt

# From my data I fitted a Von-Mises distribution, calculating Mu and Kappa. 
mu = -0.343
kappa = 10.432

# Construct random Von-Mises distribution based on Mu and Kappa values
r = np.random.vonmises(mu, kappa, 1000)

# Adjust Von-Mises curve from fitted data
x = np.linspace(-np.pi, np.pi, num=501)
y = np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa)) 

# Adjuste x limits and labels
plt.xlim(-np.pi, np.pi)
plt.xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi], 
labels=[r'$-\pi$ (0º)', r'$-\frac{\pi}{2}$ (90º)', '0 (180º)', r'$\frac{\pi}{2}$ (270º)', r'$\pi$'])

# Plot adjusted Von-Mises function as line
plt.plot(x, y, linewidth=2, color='red', zorder=3

# Plot distribution 
plt.hist(r, density=True,  bins=20, alpha=1, edgecolor='white');
plt.title('Slaty Cleavage Strike', fontweight='bold', fontsize=14)

我尝试根据问题绘制圆形直方图，例如python 中的圆形/极坐标直方图

# From the data above (mu, kappa, x and y):

theta = np.linspace(-np.pi, np.pi, num=50, endpoint=False)
radii = np.exp(kappa*np.cos(theta-mu))/(2*np.pi*i0(kappa)) 

# Bin width?
width = (2*np.pi) / 50

# Construct ax with polar projection
ax = plt.subplot(111, polar=True)

# Set Zero to North
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)

# Plot bars:
bars = ax.bar(x = theta, height = radii, width=width)
# Plot Line:
line = ax.plot(x, y, linewidth=1, color='red', zorder=3) 
 
# Grid settings
ax.set_rgrids(np.arange(1, 1.6, 0.5), angle=0, weight= 'black');

笔记：

• 我的圆形直方图在错误的方向绘制了我的数据，相差 180 度：比较两个直方图。见编辑 1
• 我相信这可能与定义在 [-pi,pi] https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.vonmises.html 上的 scipy.stats.vonmises有关编辑 1
• 我的数据最初在 [0,2pi] 之间变化。我转换为 [-pi,pi] 以拟合 Von Mises 并计算 Mu 和 Kappa。见编辑 1

我真的很想将我的数据绘制为第一个图。我的数据是地质定向数据（方位角）。有人有什么想法吗？附注。抱歉，帖子太长了。我希望这至少有帮助

编辑 1：

通过评论，我意识到有些人对数据感到困惑，无论是从[0,2pi]还是[-pi,pi]. 我意识到我的圆形直方图中绘制的错误方向来自以下内容：

1. 我的原始数据（地质数据）介于 [0,2pi]，即 0 到 360 度；
2. 然而，scipy.stats.vonmises 计算概率密度函数在 [-pi, pi]；
3. 我从我的数据中减去 pi，以便使用 scipy.stats.vonmises my_data - pi；
4. 一旦Mu和Kappa计算（CORRECTELY），我添加pi的Mu价值，以恢复原来的方向，现在之间再度范围[0,2pi]。
5. 现在，我的数据正确地指向东南：

# Add pi to fitted Mu. 
mu = - 0.343 + np.pi
kappa = 10.432

x = np.linspace(-np.pi, np.pi, num=501)
y = np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa)) 

theta = np.linspace(-np.pi, np.pi, num=50, endpoint=False)
radii = np.exp(kappa*np.cos(theta-mu))/(2*np.pi*i0(kappa)) 

# Bin width?
width = (2*np.pi) / 50

ax = plt.subplot(111, polar=True)

# Angles increase clockwise from North
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)

bars = ax.bar(x = theta, height = radii, width=width)

line = ax.plot(x, y, linewidth=1, color='red', zorder=3) 

ax.set_rgrids(np.arange(1, 1.6, 0.5), angle=0, weight= 'black');

编辑 2

正如接受的答案的评论所建议的，技巧正在改变y_lim，如下：

# SE DIRECTION
mu = - 0.343 + np.pi
kappa = 10.432


x = np.linspace(-np.pi, np.pi, num=501)
y = np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa)) 

theta = np.linspace(-np.pi, np.pi, num=50, endpoint=False)
radii = np.exp(kappa*np.cos(theta-mu))/(2*np.pi*i0(kappa)) 

# PLOT
plt.figure(figsize=(5,5))
ax = plt.subplot(111, polar=True)

# Bin width?
width = (2*np.pi) / 50

# Angles increase clockwise from North
ax.set_theta_zero_location('N'); ax.set_theta_direction(-1);

bars = ax.bar(x=theta, height = radii, width=width, bottom=0)

# Plot Line
line = ax.plot(x, y, linewidth=2, color='firebrick', zorder=3 ) 

# 'Trick': This will display Zero as a circle. Fitted Von-Mises function will lie along zero.
ax.set_ylim(-0.5, 1.5);

ax.set_rgrids(np.arange(0, 1.6, 0.5), angle=60, weight= 'bold',
             labels=np.arange(0,1.6,0.5));

最后说明：

• 提供的直方图已经标准化，因此可以用相同比例的 vM 分布绘制。

Answer 1

这就是我所取得的成就：

\n

\n

我不完全确定您是否希望 x 的范围为[-pi,pi]或[0,2pi]。如果您想要范围[0,2pi]，只需注释掉这些行ax.set_xlim并ax.set_xticks。

\n

from scipy.special import i0\nimport numpy as np\nimport matplotlib.pyplot as plt\n\n\n# From my data I fitted a Von-Mises distribution, calculating Mu and Kappa.\nmu = -0.343\nkappa = 10.432\n\n# Construct random Von-Mises distribution based on Mu and Kappa values\nr = np.random.vonmises(mu, kappa, 1000)\n\n# Adjust Von-Mises curve from fitted data\nx = np.linspace(-np.pi, np.pi, num=501)\ny = np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa))\n\n# Adjuste x limits and labels\nplt.xlim(-np.pi, np.pi)\nplt.xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi],\nlabels=[r\'$-\\pi$ (0\xc2\xba)\', r\'$-\\frac{\\pi}{2}$ (90\xc2\xba)\', \'0 (180\xc2\xba)\', r\'$\\frac{\\pi}{2}$ (270\xc2\xba)\', r\'$\\pi$\'])\n\n# Plot adjusted Von-Mises function as line\nplt.plot(x, y, linewidth=2, color=\'red\', zorder=3)\n\n# Plot distribution\nplt.hist(r, density=True,  bins=20, alpha=1, edgecolor=\'white\')\nplt.title(\'Slaty Cleavage Strike\', fontweight=\'bold\', fontsize=14)\n\n\n\n# From the data above (mu, kappa, x and y):\n\ntheta = np.linspace(-np.pi, np.pi, num=50, endpoint=False)\nradii = np.exp(kappa * np.cos(theta - mu)) / (2 * np.pi * i0(kappa))\n\n# Display width\nwidth = (2 * np.pi) / 50\n\n# Construct ax with polar projection\nax = plt.subplot(111, polar=True)\n\n# Set Orientation\nax.set_theta_zero_location(\'E\')\nax.set_theta_direction(-1)\nax.set_xlim(-np.pi/1.000001, np.pi/1.000001)  # workaround for a weird issue\nax.set_xticks([-np.pi/1.000001 + i/8 * 2*np.pi/1.000001 for i in range(8)])\n\n# Plot bars:\nbars = ax.bar(x=theta, height=radii, width=width)\n# Plot Line:\nline = ax.plot(x, y, linewidth=1, color=\'red\', zorder=3)\n\n# Grid settings\nax.set_rgrids(np.arange(.5, 1.6, 0.5), angle=0, weight=\'black\')\n\n\nplt.show()\n

\n