IMinuit¶
Minuit - программа численной минимизации функций многих переменных, широко применяемая в физике элементарных частиц. Есть два питонских интерфейса, PyMinuit и IMinuit (он особенно удобен в ipython).
from iminuit import Minuit
%pylab inline
Простой пример¶
Определим квадратичную функцию от двух параметров.
def f(a,b):
return 10*a**2+10*b**2-16*a*b+12*a-24*b
Создадим объект класса Minuit. a и b - грубые догадки, около чего надо искать минимум; error_a и error_b - оценки точности этих догадок (в начале минимизации программа будет делать шаги порядка этих величин, потом они будут уменьшаться). Пределы изменения задавать не обязательно. Валичина errordef показывает, насколько функция должна быть выше своего минимума, чтобы это считалось отклонением на одну сигму; поскольку минимизируемая функция - это, как правило, $\chi^2$, значение 1 по умолчанию вполне годится.
m=Minuit(f,a=0,error_a=1,limit_a=(-10,10),
b=0,error_b=1,limit_b=(-10,10))
Наиболее популярный метод минимизации - migrad.
m.migrad()
Значения параметров.
m.values
Значение функции в точке минимума.
m.fval
Ошибки параметров.
m.errors
Если, скажем, a - наш окончательный физический результат, то мы напишем в статье $a=1\pm0.5$. На самом деле у нас есть больше информации, поскольку ошибки a и b сильно скоррелированы. Матрица корреляции ошибок:
m.matrix()
Минимизация квадратичной формы сводится к решению системы линейных уравнений, а матрица корреляции ошибок - обратная матрица этой системы. В таком простом случае не имеет смысла использовать инструмент минимизации произвольных функций, такой, как Minuit.
M=array([[10.,-8.],[-8.,10.]])
M=inv(M)
M
M.dot(array([[-6],[12]]))
Нарисуем контуры, соответствующие отклонению на 1, 2 и 3 сигмы от оптимальной точки.
m.draw_mncontour('a','b',nsigma=3)
То же в виде цветов.
a,b,g,r=m.mncontour_grid('a','b',nsigma=3)
pcolormesh(a,b,g)
colorbar()
Дайте мне 3 параметра, и я профитирую слона. С 4 параметрами он будет махать хоботом.¶
Пусть у нас есть экспериментальные данные, и мы хотим профитировать их прямой.
def fit(a,b,x):
return a*x+b
Данные не настоящие, а сгенерированные. Все имеют ошибки 0.1.
x=linspace(0,1,11)
dy=0.1*ones(11)
y=x+dy*normal(size=11)
Функция $\chi^2$.
def chi2(a,b):
global x,y,dy
return (((y-fit(a,b,x))/dy)**2).sum()
Минимизируем.
m=Minuit(chi2,a=0,b=0,error_a=1,error_b=1)
m.migrad()
m.values
m.fval
m.matrix()
m.draw_mncontour('a','b',nsigma=3)
a,b,g,r=m.mncontour_grid('a','b',nsigma=3)
pcolormesh(a,b,g)
colorbar()
Нарисуем на одном графике экспериментальные точки, наш фит (сплошная линия) и истинную теоретическую кривую (пунктир).
errorbar(x,y,dy,fmt='ro')
xt=linspace(0,1,101)
plot(xt,fit(m.values['a'],m.values['b'],xt),'b-')
plot(xt,fit(1,0,xt),'g--')
Когда фитирующая функция есть линейная комбинация каких-то фиксированных функций с неизвестными коэффициентами, минимизация $\chi^2$ сводится к решению системы линейных уравнений. Нет надобности использовать Minuit.
Резонанс без фона¶
Пусть теперь наша фитирующая функция - Брейт-Вигнеровский резонанс (без фона), с двумя параметрами - положением и шириной (лучше бы ввести третий - высоту, но я не стал этого делать для простоты). Теперь $\chi^2$ - сложная нелинейная функция параметров.
def fit(x0,Gamma,x):
return 1/((x-x0)**2+Gamma**2)
Вот наши экспериментальные данные (с ошибками 0.1).
x=linspace(-3,3,21)
dy=0.1*ones(21)
y=fit(0,1,x)+dy*normal(size=21)
Минимизируем $\chi^2$.
def chi2(x0,Gamma):
global x,y,dy
return (((y-fit(x0,Gamma,x))/dy)**2).sum()
m=Minuit(chi2,x0=0,error_x0=1,Gamma=1,error_Gamma=1)
m.migrad()
m.values
m.fval
m.errors
m.matrix()
m.draw_mncontour('x0','Gamma',nsigma=3)
x0,Gamma,g,r=m.mncontour_grid('x0','Gamma',nsigma=3)
pcolormesh(x0,Gamma,g)
colorbar()
Теперь контуры постоянной высоты $\chi^2$ - уже не симметричные эллипсы с центром в оптимальной точке, а какие-то сложные кривые. Ошибки положения и ширины резонанса довольно-таки независимы.
Нарисуем на одном графике экспериментальные точки, наш фит (сплошная линия) и истинную теоретическую кривую (пунктир).
errorbar(x,y,dy,fmt='ro')
xt=linspace(-3.5,3.5,101)
plot(xt,fit(m.values['x0'],m.values['Gamma'],xt),'b-')
plot(xt,fit(0,1,xt),'g--')
