python: чем фиттировать распределение случайных величин

Так ты покажи, что у тебя есть, мы же даже графика не видим.

Там довольно хитрая форма, вот гистограммы: https://pix.academ.info/images/img/2017/12/11/8c3e20b73f6340abd80c8ba28574c53.... Там по физике, в принципе может быть 2 горба ну и что-то похожее видно.

Экспонента+гаусс. Никаких других распределений в физике нет. Ну разве что Пуассон и Ландау. Можешь поискать ещё Новосибирскую функцию aka логарифмический гаусс — как-то так.

Тут скорее есть ещё какой-то порог, который размазан аппаратной функцией, которая скорее всего гаусс. То есть берёшь какую-то модель и делаешь свёртку с гауссом.

И да смотри в сторону RooFit — для сложных функций лучше вычислить chi^2 и минимизировать его через RooFit.

При помощи numpy.histogram получаете координаты вершин бинов по x и y, а дальше scipy.optimize.curve_fit любой функцией на ваш вкус.

Ну я относительно недавно обсуждал подобную тему здесь.

Пробовал аппроксимировать гнуплотом (там есть стандартный fit), питоном, математикой, GSL nonlinear fit. В целом все методы работают ок, только нужно начальные значения параметров подбирать хорошо.

На картинке у тебя хороший максвелл, попробуй для начала им подобрать. С би- и мультимодальными распределениями беда, вот что нашёл.

С би- и мультимодальными распределениями беда

Да нет никакой беды.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import maxwell
from scipy.optimize import curve_fit

def maxwell_bimodal(x, n, loc1, scale1, loc2, scale2):
    P0 = maxwell.pdf(x, loc1, scale1) * (1 - n)
    P1 = maxwell.pdf(x, loc2, scale2) * n
    return(P0 + P1)

# Генерируем бимодальные данные для примера
data = np.hstack((maxwell.rvs(size=1000, loc=1, scale=1), maxwell.rvs(size=300, loc=3, scale=1))) 

binwidth = .2
bins = np.arange(0, data.max(), binwidth)
counts, binedges = np.histogram(data, bins=bins, normed=True)
bincenters = .5 * (binedges[1:] + binedges[:-1])
p0 = [.5, 1, 2, 5, .5] # Заменить на подходящие начальные параметры
popt, pcov = curve_fit(maxwell_bimodal, bincenters, counts, p0=p0)

plt.figure()
plt.hist(data, bins=bins, normed=True, color='black', alpha=.5)
plt.plot(bincenters, maxwell_bimodal(bincenters, *popt), color='red')
plt.fill(bincenters, (1 - popt[0]) * maxwell.pdf(bincenters, *popt[1:3]), alpha=.3, color='red')
plt.fill(bincenters, popt[0] * maxwell.pdf(bincenters, *popt[3:]), alpha=.3, color='red')
plt.show()

Беда не зафиттить, а подобрать подходящую функцию и объяснить, почему так, а не иначе.

В постановке задачи ТС говорит что этот этап уже пройден. А так-то, да, чем больше параметров, тем больше шанс что на данные красиво ляжет совершенно левая функция, которая к природе процесса никаким боком не относится.

А, я пропустил этот момент. Тогда проблем нет, но лучше набрать больше данных, потому что вторая мода может быть просто артефактом.

Согласен, в данном случае на это похоже. И потом, больше ещё никогда ни одной работе не навредило. ☺

Кстати, а как в таких случаях искать подходящие начальные параметры? Скажем, есть тысяча датасетов, максимумы достаточно сильно плавают, нужно всю тысячу обработать. Руками же не вариант делать.

Кстати, а как в таких случаях искать подходящие начальные параметры?

Зависит от ситуации. Например, предварительно подгонять упрощённой функцией или вычислять параметры из физической модели или из параметров самой гистограммы — число событий, RMS, максимум и т.д. Руками тоже вполне вариант — пока схема не нащупается или подгонок нужно немного.

Наиболее оптимально все это не делать вручную через параметризованные представления плотности распределения, а работать прямо с самими распределениями.

https://en.wikipedia.org/wiki/OpenBUGS ну или другие какие семплеры

Пожалуй если чисто в R, то вот это https://greta-dev.github.io/greta/