Python Scipy 中的麦克斯韦分布
Maxwellian Distribution in Python Scipy
在我感兴趣的文章中,它指出数据很好地表示为麦克斯韦分布,它还提供了平均速度 (307 km/s) 和 1 西格玛不确定性 (47 km/s) 用于分发。
使用提供的值,我尝试重新生成数据,然后使用 python scipy.stats.
使其符合麦克斯韦分布
如 中所述,scipy 中的 maxwell 函数有两个输入,1) 移动 x 变量的“loc”和 2) 对应于参数“的”a”参数a" 在麦克斯韦-玻尔兹曼方程中。
在我的例子中,我没有这些参数,因此使用 wiki page 中的均值和方差 (sigma^2) 描述,我尝试计算“a”和“loc”参数。均值和西格玛参数仅取决于“a”参数。
我遇到的第一个问题是我从 Mean (a = 192.4) 和 sigma (a = 69.8) 得到的“a”参数彼此不同。
第二个问题是我不知道如何从 Mean 和 sigma 中获得准确的 loc (shift) 值。
根据分布的形状(平均速度值落在图中,查看图 2),我试图猜测“loc”值以及从 sigma 获得的“a”值(a = 69.8),我已经生成并拟合了数据。大约看起来是正确的,但我不知道我上面提到的问题的答案,我需要一些专家的指导。感谢您的帮助。
import matplotlib.pyplot as plt
import math
from scipy.stats import norm
import random
import numpy as np
import scipy.optimize
from scipy.stats import maxwell
samplesize = 100000
mean = 307
sigma = 47
loc = 175 #my guess
a_value = np.sqrt((sigma**2 * math.pi)/(3*math.pi - 8)) #calculated based on wiki description
fig, axs = plt.subplots(1)
v_2d = maxwell.rvs(loc, a_value, size=samplesize) #array corresponding to 2D proper motion obtained from Hubbs
mean, var, skew, kurt = maxwell.stats(moments='mvsk')
N, bins, patches = plt.hist(v_2d, bins=100, density=True, alpha=0.5, histtype='bar', ec='black')
maxx = np.linspace(min(v_2d), max(v_2d), samplesize)
axs.plot(maxx, maxwell.pdf(maxx, loc, a_value), color=colorset[6], lw=2, label= r'$\mathdefault{\mu}$ = '+'{:0.1f}'.format(mean)+r' , '+r'$\mathdefault{\sigma}$ = '+'{:0.1f}'.format(sigma))
axs.set(xlabel=r'2-D Maxwellian speed (km s$^{-1}$)')
axs.set(ylabel='Frequency')
plt.legend(loc='upper right')
嗯,平均值受位置影响,而西格玛则不会。
所以从 sigma 计算 a
,计算均值,就好像 loc=0,找到差异并将其分配给位置,对 100K RV 进行采样以检查是否
采样 mean/stddev 足够接近。
代码,Python3.8,Windows10 x64
import numpy as np
from scipy.stats import maxwell
σ = 47
μ = 307
a = σ * np.sqrt(np.pi/(3.0*np.pi - 8.0))
print(a)
m = 2.0*a*np.sqrt(2.0/np.pi)
print(m) # as if loc=0
loc = μ - m
print(loc)
print("----------Now test--------------------")
# sampling
q = maxwell.rvs(loc=loc, scale=a, size=100000)
print(np.mean(q))
print(np.std(q))
作为我的输出
306.9022249667151
47.05319429681308
够好吗?
在我感兴趣的文章中,它指出数据很好地表示为麦克斯韦分布,它还提供了平均速度 (307 km/s) 和 1 西格玛不确定性 (47 km/s) 用于分发。
使用提供的值,我尝试重新生成数据,然后使用 python scipy.stats.
使其符合麦克斯韦分布如
在我的例子中,我没有这些参数,因此使用 wiki page 中的均值和方差 (sigma^2) 描述,我尝试计算“a”和“loc”参数。均值和西格玛参数仅取决于“a”参数。
我遇到的第一个问题是我从 Mean (a = 192.4) 和 sigma (a = 69.8) 得到的“a”参数彼此不同。 第二个问题是我不知道如何从 Mean 和 sigma 中获得准确的 loc (shift) 值。
根据分布的形状(平均速度值落在图中,查看图 2),我试图猜测“loc”值以及从 sigma 获得的“a”值(a = 69.8),我已经生成并拟合了数据。大约看起来是正确的,但我不知道我上面提到的问题的答案,我需要一些专家的指导。感谢您的帮助。
import matplotlib.pyplot as plt
import math
from scipy.stats import norm
import random
import numpy as np
import scipy.optimize
from scipy.stats import maxwell
samplesize = 100000
mean = 307
sigma = 47
loc = 175 #my guess
a_value = np.sqrt((sigma**2 * math.pi)/(3*math.pi - 8)) #calculated based on wiki description
fig, axs = plt.subplots(1)
v_2d = maxwell.rvs(loc, a_value, size=samplesize) #array corresponding to 2D proper motion obtained from Hubbs
mean, var, skew, kurt = maxwell.stats(moments='mvsk')
N, bins, patches = plt.hist(v_2d, bins=100, density=True, alpha=0.5, histtype='bar', ec='black')
maxx = np.linspace(min(v_2d), max(v_2d), samplesize)
axs.plot(maxx, maxwell.pdf(maxx, loc, a_value), color=colorset[6], lw=2, label= r'$\mathdefault{\mu}$ = '+'{:0.1f}'.format(mean)+r' , '+r'$\mathdefault{\sigma}$ = '+'{:0.1f}'.format(sigma))
axs.set(xlabel=r'2-D Maxwellian speed (km s$^{-1}$)')
axs.set(ylabel='Frequency')
plt.legend(loc='upper right')
嗯,平均值受位置影响,而西格玛则不会。
所以从 sigma 计算 a
,计算均值,就好像 loc=0,找到差异并将其分配给位置,对 100K RV 进行采样以检查是否
采样 mean/stddev 足够接近。
代码,Python3.8,Windows10 x64
import numpy as np
from scipy.stats import maxwell
σ = 47
μ = 307
a = σ * np.sqrt(np.pi/(3.0*np.pi - 8.0))
print(a)
m = 2.0*a*np.sqrt(2.0/np.pi)
print(m) # as if loc=0
loc = μ - m
print(loc)
print("----------Now test--------------------")
# sampling
q = maxwell.rvs(loc=loc, scale=a, size=100000)
print(np.mean(q))
print(np.std(q))
作为我的输出
306.9022249667151
47.05319429681308
够好吗?