Распределение Максвелла в Python Scipy - PullRequest
2 голосов
/ 07 августа 2020

В интересующей меня статье говорится, что данные хорошо представлены с распределением Максвелла, а также приводится средняя скорость (307 км / с) и неопределенность 1 сигма (47 км / с) для распределения. .

Используя предоставленные значения, я попытался повторно сгенерировать данные, а затем подогнать их к распределению Максвелла, используя python scipy.stats.

Как описано в здесь , функция maxwell в scipy принимает два входа: 1) «lo c», который сдвигает переменную x, и 2) параметр «a», который соответствует параметру «a» в уравнении Максвелла-Больцмана.

В моем случае у меня нет ни одного из этих параметров, поэтому, используя описание Среднее и дисперсия (сигма ^ 2) на вики-странице , я попытался вычислить «а» и «lo». c "параметр. Оба параметра mean и sigma зависят только от параметра «a».

Первой проблемой, с которой я столкнулся, был параметр «a», который я получил от Mean (a = 192,4) и sigma (a = 69,8). отличаются друг от друга. Вторая проблема заключается в том, что я не знаю, как получить точное значение lo c (сдвиг) из среднего и сигма.

На основе формы распределения (где значения средней скорости попадают в график, проверьте рисунок 2), я попытался угадать значение «lo c» и вместе со значением «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')

Distribution based on my guess

Средняя и медианная скорость в распределении Максвелла

1 Ответ

1 голос
/ 08 августа 2020

Ну, среднее значение зависит от местоположения, а сигма - нет. Итак, вычислите a из сигмы, вычислите среднее значение, как если бы loc = 0, найдите разницу и назначьте ее местоположению, выполните выборку 100K RV, чтобы проверить, достаточно ли близки выборочное среднее / stddev.

Код, Python 3.8, Windows 10 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

Достаточно хорошо?

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...