У меня есть два массива с данными x и y.
Эти данные показывают ненормальное поведение. Мне нужен график соответствия, а также мю и сигма, чтобы сделать некоторые статистические данные.
Я подобрал, чтобы вычислить мю, сигму и, далее, некоторые ее статистические значения. (См. Код ниже)
Я получаю масштабный коэффициент, с помощью которого мне нужно умножить распределение на интеграл по точкам данных.
Код ниже, работает. Мой вопрос сейчас заключается в том, если (я уверен) есть лучший способ сделать это? Это похоже на обходной путь, который иногда будет работать. Я хочу лучший способ сделать это, потому что мне нужно подготовить сотни таких.
Мой код (извините, что так долго, хотел включить все, кроме импорта сырых данных):
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# produce plot True/False
ploton = True
x0=np.array([3.58381e+01, 3.27125e+01, 2.98680e+01, 2.72888e+01, 2.49364e+01,
2.27933e+01, 2.08366e+01, 1.90563e+01, 1.74380e+01, 1.59550e+01,
1.45904e+01, 1.33460e+01, 1.22096e+01, 1.11733e+01, 1.02262e+01,
9.35893e+00, 8.56556e+00, 7.86688e+00, 7.20265e+00, 6.59782e+00,
6.01571e+00, 5.53207e+00, 5.03979e+00, 4.64415e+00, 4.19920e+00,
3.83595e+00, 3.50393e+00, 3.28070e+00, 3.00930e+00, 2.75634e+00,
2.52050e+00, 2.31349e+00, 2.12280e+00, 1.92642e+00, 1.77820e+00,
1.61692e+00, 1.49094e+00, 1.36233e+00, 1.22935e+00, 1.14177e+00,
1.03078e+00, 9.39603e-01, 8.78425e-01, 1.01490e+00, 1.07461e-01,
4.81523e-02, 4.81523e-02, 1.00000e-02, 1.00000e-02])
y0=np.array([3.94604811e+04, 2.78223936e+04, 1.95979179e+04, 2.14447807e+04,
1.68677487e+04, 1.79429516e+04, 1.73589776e+04, 2.16101026e+04,
3.79705638e+04, 6.83622301e+04, 1.73687772e+05, 5.74854475e+05,
1.69497465e+06, 3.79135941e+06, 7.76757753e+06, 1.33429094e+07,
1.96096415e+07, 2.50403065e+07, 2.72818618e+07, 2.53120387e+07,
1.93102362e+07, 1.22219224e+07, 4.96725699e+06, 1.61174658e+06,
3.19352386e+05, 1.80305856e+05, 1.41728002e+05, 1.66191809e+05,
1.33223816e+05, 1.31384905e+05, 2.49100945e+05, 2.28300583e+05,
3.01063903e+05, 1.84271914e+05, 1.26412781e+05, 8.57488083e+04,
1.35536571e+05, 4.50076293e+04, 1.98080100e+05, 2.27630303e+05,
1.89484527e+05, 0.00000000e+00, 1.36543525e+05, 2.20677520e+05,
3.60100586e+05, 1.62676486e+05, 1.90105093e+04, 9.27461467e+05,
1.58373542e+05])
Dnm = x0
dndlndp = y0
#lognormal PDF:
def f(x, mu, sigma) :
return 1/(np.sqrt(2*np.pi)*sigma*x)*np.exp(-((np.log(x)-mu)**2)/(2*sigma**2))
#normalizing y-values to obtain lognormal distributed data:
y0_normalized = y0/np.trapz(x0.ravel(), y0.ravel())
#calculating mu/sigma of this distribution:
params, extras = curve_fit(f, x0.ravel(), y0_normalized.ravel())
median = np.exp(params[0])
mu = params[0]
sigma = params[1]
#output of mu / sigma / calculated median:
print "mu=%g, sigma=%g" % (params[0], params[1])
print "median=%g" % median
#new variable z for smooth fit-curve:
z = np.linspace(0.1, 100, 10000)
#######################
Dnm = np.ravel(Dnm)
dndlndp = np.ravel(dndlndp)
Dnm_rev = list(reversed(Dnm))
dndlndp_rev = list(reversed(dndlndp))
scalingfactor = np.trapz(dndlndp_rev, Dnm_rev, dx = np.log(Dnm_rev))
#####################
#plotting
if ploton:
plt.plot(z, f(z, mu, sigma)*scalingfactor, label="fit", color = "red")
plt.scatter(x0, y0, label="data")
plt.xlim(3,20)
plt.xscale("log")
plt.legend()
РЕДАКТИРОВАТЬ1: Может быть, я должен добавить, что я понятия не имею, почему коэффициент масштабирования рассчитывается с
scalingfactor = np.trapz(dndlndp_rev, Dnm_rev, dx = np.log(Dnm_rev))
верно. Это была просто попытка и ошибка. Я действительно хочу знать, почему это делает трюк, так как "область" всех объединенных бункеров:
N = np.trapz(dndlndp_rev, np.log(Dnm_rev), dx = np.log(Dnm_rev))
потому что ширина ящиков равна log (Dnm).
EDIT2: Спасибо за все ответы. Я скопировал массивы в код, который теперь можно запустить. Я хочу упростить вопрос, так как думаю, что из-за моего плохого английского я не смог сказать, чего я действительно хочу:
У меня есть ненормальный набор данных. Код выше позволяет мне рассчитать мю и сигму. Для этого мне нужно нормализовать данные, и область под функцией теперь = 1.
Чтобы построить логнормальную функцию с вычисленными значениями mu и sigma, мне нужно умножить функцию на (неизвестный) коэффициент, потому что область под реальной функцией - что-то вроде 1e8, но не одно. Я сделал обходной путь, вычислив этот «коэффициент масштабирования» через интеграл trapz дискретных необработанных данных.
Должен быть лучший способ для построения подобранной функции, когда му и сигма уже известны.