mises分布的python圆直方图

mwngjboj  于 2021-07-13  发布在  Java
关注(0)|答案(1)|浏览(355)

在过去的几天里,我一直在尝试用python绘制循环数据,通过构建一个从0到2pi的循环直方图并拟合von mises分布。我真正想要实现的是:
拟合von mises分布的定向数据。此绘图由matplotlib、scipy和numpy构建,可在以下位置找到:http://jpktd.blogspot.com/2012/11/polar-histogram.html

这个图是用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
我相信这可能与scipy.stats.vonmises有关,它是在[-pi,pi]上定义的https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.vonmises.htmlsee 编辑1
我的数据最初在[0,2pi]之间变化。我转换成[-pi,pi]来拟合von mises并计算mu和kappa。请参见编辑1
我真的很想把我的数据作为第一个绘图之一。我的数据是地质方向数据(方位角)。有人有什么想法吗?另外,很抱歉发了这么长的邮件。我希望这至少有帮助

编辑1:

通过浏览这些评论,我意识到有一些人对数据是否从 [0,2pi] 或者 [-pi,pi] . 我意识到在圆形直方图中绘制的错误方向来自以下几点:
我的原始数据(地质数据)介于 [0,2pi] ,即0至360度;
但是,scipy.stats.vonmises在 [-pi, pi] ;
我从数据中减去pi,以便使用scipy.stats.vonmises my_data - pi ;
一次 Mu 以及 Kappa 我补充说,这是经过计算的(正确的) piMu 值,恢复原来的方向,现在再次介于 [0,2pi] .
现在,我的数据正确地指向东南:


# 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));

要做的工作:

要么规范化直方图,要么将von mises分布乘以样本的数量,如注解所示。

vwkv1x7d

vwkv1x7d1#

这就是我的成就:

我不太确定你是否希望x的范围是 [-pi,pi] 或者 [0,2pi] . 如果你想要这个范围 [0,2pi] 相反,只需注解掉这些行 ax.set_xlim 以及 ax.set_xticks .

from scipy.special import i0
import numpy as np
import matplotlib.pyplot 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)

# 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))

# Display width

width = (2 * np.pi) / 50

# Construct ax with polar projection

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

# Set Orientation

ax.set_theta_zero_location('E')
ax.set_theta_direction(-1)
ax.set_xlim(-np.pi/1.000001, np.pi/1.000001)  # workaround for a weird issue
ax.set_xticks([-np.pi/1.000001 + i/8 * 2*np.pi/1.000001 for i in range(8)])

# 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(.5, 1.6, 0.5), angle=0, weight='black')

plt.show()

相关问题