using Pkg
Pkg.add(url="https://github.com/fkguo/IMinuit.jl")
Updating `~/IMinuit.jl/docs/Project.toml` [beb75e20] - IMinuit v0.2.0 `~/IMinuit.jl` Updating `~/IMinuit.jl/docs/Manifest.toml` [861a8166] - Combinatorics v1.0.2 [bbf7d656] - CommonSubexpressions v0.3.0 [163ba53b] - DiffResults v1.1.0 [b552c78f] - DiffRules v1.15.1 [f6369f11] - ForwardDiff v0.10.35 [beb75e20] - IMinuit v0.2.0 `~/IMinuit.jl` [1e83bf80] - StaticArraysCore v1.4.2 [8ba89e20] - Distributed Resolving package versions... Updating `~/IMinuit.jl/docs/Project.toml` [beb75e20] + IMinuit v0.2.0 `~/IMinuit.jl` Updating `~/IMinuit.jl/docs/Manifest.toml` [861a8166] + Combinatorics v1.0.2 [bbf7d656] + CommonSubexpressions v0.3.0 [163ba53b] + DiffResults v1.1.0 [b552c78f] + DiffRules v1.15.1 [f6369f11] + ForwardDiff v0.10.35 [beb75e20] + IMinuit v0.2.0 `~/IMinuit.jl` [1e83bf80] + StaticArraysCore v1.4.2 [8ba89e20] + Distributed
using IMinuit
using BenchmarkTools
[ Info: Precompiling IMinuit [beb75e20-2205-47e6-ad51-640e9c2309f1]
initiated!!
f(x) = x[1]^2 + (x[2]-1)^2 + (x[3]-2)^4
f1(x, y, z) = x^2 + (y-1)^2 + (z-2)^4
f1 (generic function with 1 method)
# using array parameters
m = Minuit(f, [1, 1, 4], fix_x0 = true)
migrad(m)
Migrad | ||||
---|---|---|---|---|
FCN = 1 | Nfcn = 45 | |||
EDM = 2.47e-05 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | x0 | 1.00 | 0.01 | yes | ||||
1 | x1 | 1 | 1 | |||||
2 | x2 | 2 | 5 |
x0 | x1 | x2 | |
---|---|---|---|
x0 | 0 | 0 | 0 |
x1 | 0 | 1 | -1.22e-09 |
x2 | 0 | -1.22e-09 | 27.4 |
# a new fit can continue from the previous fit
m_new = Minuit(f, m, fix_x0 = false)
migrad(m_new)
Migrad | ||||
---|---|---|---|---|
FCN = 7.285e-06 | Nfcn = 40 | |||
EDM = 4.85e-06 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | x0 | -0 | 1 | |||||
1 | x1 | 1 | 1 | |||||
2 | x2 | 2 | 8 |
x0 | x1 | x2 | |
---|---|---|---|
x0 | 1 | 1.78e-15 | 4.96e-29 |
x1 | 1.78e-15 | 1 | 2.79e-14 |
x2 | 4.96e-29 | 2.79e-14 | 61.7 |
# using array parameters, using `ForwardDiff: gradient` to compute the gradient
gradf(x) = gradient(f, x)
mgrad = Minuit(f, [1, 1, 4], grad = gradf)
migrad(mgrad)
Migrad | ||||
---|---|---|---|---|
FCN = 0.0001683 | Nfcn = 59, Ngrad = 6 | |||
EDM = 0.000115 (Goal: 0.0002) | time = 0.6 sec | |||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | x0 | 0 | 1 | |||||
1 | x1 | 1 | 1 | |||||
2 | x2 | 2 | 4 |
x0 | x1 | x2 | |
---|---|---|---|
x0 | 1 | 0 | 0 |
x1 | 0 | 1 | 0 |
x2 | 0 | 0 | 13.2 |
# parameters are given individually
m1 = Minuit(f1, x = 1, y = 1, z = 4)
migrad(m1)
Migrad | ||||
---|---|---|---|---|
FCN = 0.0001683 | Nfcn = 89 | |||
EDM = 0.000115 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | x | 0 | 1 | |||||
1 | y | 1 | 1 | |||||
2 | z | 2 | 4 |
x | y | z | |
---|---|---|---|
x | 1 | 0 | 0 |
y | 0 | 1 | 0 |
z | 0 | 0 | 13.2 |
@show iminuit.__version__
propertynames(m)
iminuit.__version__ = "2.18.0"
101-element Vector{Symbol}: :LEAST_SQUARES :LIKELIHOOD :__annotations__ :__class__ :__delattr__ :__dir__ :__doc__ :__eq__ :__format__ :__ge__ :__getattribute__ :__gt__ :__hash__ ⋮ :profile :reset :scan :scipy :simplex :strategy :throw_nan :tol :valid :values :var2pos :visualize
# the doc strings are from `iminuit`
@doc migrad
Docstring pulled from the Python `iminuit`: Run Migrad minimization. Migrad from the Minuit2 library is a robust minimisation algorithm which earned its reputation in 40+ years of almost exclusive usage in high-energy physics. How Migrad works is described in the `Minuit paper`_. It uses first and approximate second derivatives to achieve quadratic convergence near the minimum. Parameters ---------- ncall : Approximate maximum number of calls before minimization will be aborted. If set to None, use the adaptive heuristic from the Minuit2 library (Default: None). Note: The limit may be slightly violated, because the condition is checked only after a full iteration of the algorithm, which usually performs several function calls. iterate : Automatically call Migrad up to N times if convergence was not reached (Default: 5). This simple heuristic makes Migrad converge more often even if the numerical precision of the cost function is low. Setting this to 1 disables the feature. See Also -------- simplex, scan
The data are taken from BES Collaboration, Phys. Rev. D 62 (2000) 032002.
Here we use a simple model, which is not meant to be the correct one, to fit to the data.
using CSV
using DataFrames
using Plots
pyplot(framestyle = :box, minorticks = 5)
using LaTeXStrings
[ Info: Precompiling CSV [336ed68f-0bac-5ca0-87d4-7b16caf5d00b] [ Info: Precompiling IJuliaExt [2f4121a4-3b3a-5ce6-9c5e-1f2673ce168a]
data_df = DataFrame(CSV.File("./testdata.csv"))
const data = Data(data_df)
Data([0.303, 0.309, 0.321, 0.327, 0.333, 0.339, 0.345, 0.351, 0.357, 0.363 … 0.531, 0.537, 0.543, 0.549, 0.555, 0.561, 0.567, 0.573, 0.579, 0.585], [8.522656, 44.87459, 29.63286, 58.13258, 28.2143, 129.0572, 181.7821, 199.5301, 269.4043, 186.1139 … 5217.45, 5867.647, 5717.979, 5527.38, 5548.063, 5386.82, 5425.564, 4744.075, 3899.626, 2725.864], [14.44523, 78.24785, 23.01279, 28.29703, 14.79316, 35.48893, 43.7948, 42.20547, 49.41916, 42.54225 … 164.0276, 174.4779, 169.0323, 165.0818, 164.632, 160.1847, 160.5301, 149.5817, 138.9569, 117.8775], 47)
@plt_data data xlab=L"m_{\pi\pi}"*" [GeV]" ylab="Events"
sys:1: UserWarning: You passed a edgecolor/edgecolors ((0.0, 0.0, 0.0, 1.0)) for an unfilled marker ('_'). Matplotlib is ignoring the edgecolor in favor of the facecolor. This behavior may change in the future. sys:1: UserWarning: No data for colormapping provided via 'c'. Parameters 'vmin', 'vmax' will be ignored
const M = 3.686; const mπ = 0.14; const mJ = 3.097;
λ(x, y, z) = x^2 + y^2 + z^2 - 2x*y - 2y*z - 2z*x
# a simple function that will be used to fit the data: QCD multipole expansion model for ψ'→J/ψπ⁺π⁻
# The important ππ FSI effect is not taken into account
# bg is just for introducing a third parameter
function dist(w, N, c, bg)
if (w ≤ 2mπ || w ≥ M-mJ)
res = 0.0
else
q1 = sqrt(λ(w^2, mπ^2, mπ^2))/(2w)
q2 = sqrt(λ(M^2, w^2, mJ^2))/(2M)
res = N * q1 * q2 * (w^2 - c*mπ^2)^2 + bg
end
return res * 1e6
end;
dist(x, p) = dist(x, p...)
dist (generic function with 2 methods)
# parameters given individually
χsq1(N, c, bg) = chisq(dist, data, (N, c, bg));
# all parameters are vairables of χsq
fit1 = Minuit(χsq1, N = 1, c = 2, bg = 0, error_N = 0.1, error_c = 0.1, error_bg = 0.1)
fit1.strategy = 1;
# parameters are collected into a tuple or an array, which is the only variable of χsq
parname = [:N, :c, :bg]
χsq(par) = chisq(dist, data, par)
gradf(par) = gradient(χsq, par)
fit = Minuit(χsq, [1, 2, 0], error = 0.1*ones(3), name = parname, grad = gradf)
fit.strategy = 1;
# or simply using model_fit or @model_fit
fit2 = model_fit(dist, data, [1, 2, 0], error = 0.1*ones(3), name = parname, fix_bg = true)
fit2.strategy = 1;
migrad(fit2)
Migrad | ||||
---|---|---|---|---|
FCN = 68.22 | Nfcn = 60 | |||
EDM = 6.15e-05 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | N | 2.67 | 0.04 | |||||
1 | c | 4.33 | 0.06 | |||||
2 | bg | 0.0 | 0.1 | yes |
N | c | bg | |
---|---|---|---|
N | 0.00178 | 0.00233 (0.893) | 0 |
c | 0.00233 (0.893) | 0.0038 | 0 |
bg | 0 | 0 | 0 |
# the privous fit status can be passed to a new fit
fit2_new = model_fit(dist, data, fit2, name = parname, fix_bg = false)
migrad(fit2_new)
Migrad | ||||
---|---|---|---|---|
FCN = 64.68 | Nfcn = 65 | |||
EDM = 2.08e-05 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | N | 2.61 | 0.05 | |||||
1 | c | 4.20 | 0.09 | |||||
2 | bg | -0.022e-3 | 0.012e-3 |
N | c | bg | |
---|---|---|---|
N | 0.00268 | 0.00446 (0.925) | 3.73e-07 (0.605) |
c | 0.00446 (0.925) | 0.00867 | 8.29e-07 (0.748) |
bg | 3.73e-07 (0.605) | 8.29e-07 (0.748) | 1.41e-10 |
@btime migrad(fit)
minos(fit)
migrad(fit)
215.846 μs (4818 allocations: 93.65 KiB)
Migrad | ||||
---|---|---|---|---|
FCN = 64.68 | Nfcn = 399650, Ngrad = 11 | |||
EDM = 2.34e-15 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | N | 2.61 | 0.05 | -0.05 | 0.05 | |||
1 | c | 4.20 | 0.09 | -0.09 | 0.09 | |||
2 | bg | -0.022e-3 | 0.012e-3 | -0.012e-3 | 0.012e-3 |
N | c | bg | ||||
---|---|---|---|---|---|---|
Error | -0.05 | 0.05 | -0.09 | 0.09 | -0.012e-3 | 0.012e-3 |
Valid | True | True | True | True | True | True |
At Limit | False | False | False | False | False | False |
Max FCN | False | False | False | False | False | False |
New Min | False | False | False | False | False | False |
N | c | bg | |
---|---|---|---|
N | 0.00268 | 0.00446 (0.925) | 3.73e-07 (0.605) |
c | 0.00446 (0.925) | 0.00867 | 8.29e-07 (0.748) |
bg | 3.73e-07 (0.605) | 8.29e-07 (0.748) | 1.42e-10 |
@btime migrad(fit1)
minos(fit1)
migrad(fit1)
136.016 μs (730 allocations: 19.93 KiB)
Migrad | ||||
---|---|---|---|---|
FCN = 64.68 | Nfcn = 565652 | |||
EDM = 1.52e-14 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | N | 2.61 | 0.05 | -0.05 | 0.05 | |||
1 | c | 4.20 | 0.09 | -0.09 | 0.09 | |||
2 | bg | -0.022e-3 | 0.012e-3 | -0.012e-3 | 0.012e-3 |
N | c | bg | ||||
---|---|---|---|---|---|---|
Error | -0.05 | 0.05 | -0.09 | 0.09 | -0.012e-3 | 0.012e-3 |
Valid | True | True | True | True | True | True |
At Limit | False | False | False | False | False | False |
Max FCN | False | False | False | False | False | False |
New Min | False | False | False | False | False | False |
N | c | bg | |
---|---|---|---|
N | 0.00268 | 0.00446 (0.925) | 3.73e-07 (0.606) |
c | 0.00446 (0.925) | 0.00867 | 8.29e-07 (0.748) |
bg | 3.73e-07 (0.606) | 8.29e-07 (0.748) | 1.42e-10 |
minos(fit1)
Migrad | ||||
---|---|---|---|---|
FCN = 64.68 | Nfcn = 565820 | |||
EDM = 1.52e-14 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | N | 2.61 | 0.05 | -0.05 | 0.05 | |||
1 | c | 4.20 | 0.09 | -0.09 | 0.09 | |||
2 | bg | -0.022e-3 | 0.012e-3 | -0.012e-3 | 0.012e-3 |
N | c | bg | ||||
---|---|---|---|---|---|---|
Error | -0.05 | 0.05 | -0.09 | 0.09 | -0.012e-3 | 0.012e-3 |
Valid | True | True | True | True | True | True |
At Limit | False | False | False | False | False | False |
Max FCN | False | False | False | False | False | False |
New Min | False | False | False | False | False | False |
N | c | bg | |
---|---|---|---|
N | 0.00268 | 0.00446 (0.925) | 3.73e-07 (0.606) |
c | 0.00446 (0.925) | 0.00867 | 8.29e-07 (0.748) |
bg | 3.73e-07 (0.606) | 8.29e-07 (0.748) | 1.42e-10 |
# the ordering of dist, fit and data does not matter
@plt_best dist fit data
# MIGRAD contour of two parameters with the other ones fixed
# needs PyPlot or the pyplot backend of Plots
draw_contour(fit1,:N, :c, bound=3, bins=100)
([2.4561531846064364, 2.4592911671637516, 2.4624291497210673, 2.4655671322783825, 2.4687051148356978, 2.4718430973930134, 2.4749810799503287, 2.4781190625076444, 2.4812570450649596, 2.484395027622275 … 2.738571614764822, 2.7417095973221373, 2.7448475798794525, 2.747985562436768, 2.7511235449940834, 2.7542615275513986, 2.7573995101087143, 2.7605374926660295, 2.7636754752233452, 2.7668134577806605], [3.9216807759073293, 3.9273242217384023, 3.932967667569476, 3.938611113400549, 3.944254559231622, 3.9498980050626953, 3.9555414508937683, 3.9611848967248418, 3.9668283425559148, 3.972471788386988 … 4.429590900703915, 4.435234346534989, 4.440877792366061, 4.446521238197135, 4.452164684028208, 4.457808129859281, 4.463451575690354, 4.469095021521428, 4.4747384673525, 4.480381913183574], [20.970459755972996 20.02080897997878 … 290.2779850474049 296.35529928065046; 21.146256518693036 20.129488037901638 … 285.04032680263396 291.0733880873432; … ; 350.8704215520153 342.54385622437354 … 19.824213491646248 20.933801142223516; 357.4930611817346 349.0828960981299 … 19.59349826608762 20.64582394490988])
# contour of parameter space from MINOS
draw_mncontour(fit1,:N, :c, nsigma=3, numpoints=100)
PyObject <matplotlib.contour.ContourSet object at 0x7f14bf5f5c60>
matrix(fit1, correlation = true)
3×3 Matrix{Float64}: 1.0 0.925213 0.605526 0.925213 1.0 0.748397 0.605526 0.748397 1.0
@show fit1.matrix
matrix(fit1)
fit1.matrix = [0.0026808279257967637 0.004460717050849894 3.729529830296808e-7; 0.004460717050849894 0.008670748910938015 8.289868474538332e-7; 3.729529830296808e-7 8.289868474538332e-7 1.4150579833676482e-10]
3×3 Matrix{Float64}: 0.00268083 0.00446072 3.72953e-7 0.00446072 0.00867075 8.28987e-7 3.72953e-7 8.28987e-7 1.41506e-10
# this gives parameter sets at the 1σ boundary
@time contour_df(fit1, χsq1, npts = 5)
1.702498 seconds (3.74 M allocations: 246.451 MiB, 5.41% gc time, 97.76% compilation time)
Row | chisq | N | c | bg |
---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | |
1 | 64.6799 | 2.61148 | 4.20103 | -2.19212e-5 |
2 | 65.6799 | 2.56003 | 4.11414 | -2.93031e-5 |
3 | 65.6799 | 2.66359 | 4.28643 | -1.48855e-5 |
4 | 65.68 | 2.56385 | 4.10711 | -3.10184e-5 |
5 | 65.6799 | 2.65961 | 4.29336 | -1.32339e-5 |
6 | 65.68 | 2.5798 | 4.13027 | -3.39527e-5 |
7 | 65.68 | 2.64221 | 4.2688 | -1.01556e-5 |
8 | 65.7184 | 2.56003 | 4.10711 | -3.0459e-5 |
9 | 65.7184 | 2.56003 | 4.10711 | -3.0459e-5 |
10 | 65.7193 | 2.66359 | 4.29336 | -1.38007e-5 |
11 | 65.7193 | 2.66359 | 4.29336 | -1.38007e-5 |
# random sampling of parameters in given ranges, keeping those within 1σ
@time parsam_df = contour_df_samples(fit1, χsq1, (:N, :c), ([2.5,2.8], [4.0,4.3]), nsamples = 3000)
2.307346 seconds (6.62 M allocations: 264.082 MiB, 5.05% gc time, 27.12% compilation time)
Row | chisq | N | c | bg |
---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | |
1 | 64.6799 | 2.61148 | 4.20103 | -2.19212e-5 |
2 | 65.3583 | 2.65105 | 4.27719 | -1.47288e-5 |
3 | 64.7453 | 2.62234 | 4.21387 | -2.12807e-5 |
4 | 64.8934 | 2.59253 | 4.15935 | -2.62679e-5 |
5 | 64.8994 | 2.60053 | 4.19757 | -2.09793e-5 |
6 | 64.8844 | 2.59933 | 4.16706 | -2.5927e-5 |
7 | 65.2471 | 2.64545 | 4.24388 | -1.94565e-5 |
8 | 65.542 | 2.58823 | 4.19096 | -2.03662e-5 |
9 | 64.7453 | 2.62474 | 4.22317 | -2.00633e-5 |
10 | 65.3307 | 2.64605 | 4.27399 | -1.45818e-5 |
11 | 64.9274 | 2.60544 | 4.20807 | -1.9912e-5 |
12 | 65.1491 | 2.62244 | 4.24218 | -1.66078e-5 |
13 | 65.0551 | 2.62124 | 4.23778 | -1.71687e-5 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
185 | 65.0177 | 2.63875 | 4.23728 | -1.96385e-5 |
186 | 65.3285 | 2.57082 | 4.12684 | -2.86867e-5 |
187 | 65.1625 | 2.59883 | 4.20297 | -1.98453e-5 |
188 | 64.877 | 2.62464 | 4.23568 | -1.79791e-5 |
189 | 65.4006 | 2.64295 | 4.27379 | -1.41986e-5 |
190 | 64.9888 | 2.62404 | 4.23948 | -1.72703e-5 |
191 | 65.1569 | 2.57603 | 4.14365 | -2.65807e-5 |
192 | 65.0229 | 2.61664 | 4.22998 | -1.78266e-5 |
193 | 64.8011 | 2.62954 | 4.23098 | -1.94266e-5 |
194 | 65.0025 | 2.62404 | 4.20377 | -2.32039e-5 |
195 | 64.9782 | 2.58453 | 4.15015 | -2.669e-5 |
196 | 65.6403 | 2.57332 | 4.15975 | -2.34856e-5 |
# get parameter ranges
extrema(parsam_df.:N), extrema(parsam_df.:c)
((2.566922307435812, 2.660053351117039), (4.115238412804268, 4.2888962987662556))
scatter(parsam_df.:N, parsam_df.:c, xlab = "N", ylab = "c")
@time contour_df_samples(fit, χsq, (:x0, :x1, :x2), ([2.5,2.8], [4.0,4.3], (-3e-5,3e-5)), nsamples = 1000)
0.983029 seconds (1.82 M allocations: 98.122 MiB, 4.15% gc time, 74.24% compilation time)
Row | chisq | N | c | bg |
---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | |
1 | 64.6799 | 2.61148 | 4.20103 | -2.19212e-5 |
2 | 64.9362 | 2.60601 | 4.20901 | -2.06306e-5 |
3 | 65.671 | 2.62162 | 4.2006 | -1.71471e-5 |
4 | 65.137 | 2.5979 | 4.19099 | -1.79279e-5 |
5 | 65.094 | 2.6006 | 4.19159 | -1.78679e-5 |
6 | 65.3406 | 2.63333 | 4.26126 | -1.40841e-5 |
7 | 65.1599 | 2.5988 | 4.15856 | -2.54955e-5 |
@time contour_df_samples(fit, χsq, :x0, (2.5,2.8), nsamples = 20)
0.153705 seconds (550.96 k allocations: 18.102 MiB, 87.97% compilation time)
Row | chisq | N | c | bg |
---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | |
1 | 64.6799 | 2.61148 | 4.20103 | -2.19212e-5 |
2 | 65.5614 | 2.56316 | 4.11968 | -2.88302e-5 |
3 | 65.0779 | 2.57895 | 4.14646 | -2.65198e-5 |
4 | 64.785 | 2.59474 | 4.17304 | -2.42721e-5 |
5 | 64.6803 | 2.61053 | 4.19944 | -2.20544e-5 |
6 | 64.7617 | 2.62632 | 4.2256 | -1.98768e-5 |
7 | 65.0271 | 2.64211 | 4.2515 | -1.77407e-5 |
8 | 65.4744 | 2.65789 | 4.27721 | -1.56345e-5 |
The example is Example A: Fit of a gaussian model to a histogram
using PyCall
# import numpy from Python to generate the same data as in the example
np = pyimport(:numpy)
default_rng = pyimport("numpy.random").default_rng
rng = default_rng(seed=1)
const w, xe = np.histogram(rng.normal(0, 1, 10000), bins=1000)
([1, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], [-3.8378621427178974, -3.830092502694802, -3.8223228626717063, -3.8145532226486107, -3.806783582625515, -3.7990139426024196, -3.791244302579324, -3.7834746625562286, -3.775705022533133, -3.7679353825100375 … 3.8618511201697956, 3.869620760192891, 3.8773904002159867, 3.8851600402390822, 3.892929680262178, 3.9006993202852733, 3.908468960308369, 3.9162386003314644, 3.92400824035456, 3.931777880377655])
# define the model and the score function to minimize
using SpecialFunctions
function cdf(x, par)
mu, sigma = par
z = (x - mu) / sigma
return 0.5 * (1 + erf(z / sqrt(2)))
end
function score(par)
amp = par[1]
rest = par[2:end]
mu = amp * (cdf.(xe[2:end], Ref(rest)) - cdf.(xe[1:end-1], Ref(rest)) )
return 2 * sum(@. mu - w * log(mu + 1e-100))
end
score (generic function with 1 method)
const start_values = [1.5 * sum(w), 1.0, 2.0]
const limits = [(0, nothing), nothing, (0, nothing)];
# w/o grad
m = Minuit(score, start_values, limit=limits)
m.strategy = 0
# using grad
grad_fd(pars) = gradient(score, pars)
m_fd = Minuit(score, start_values, limit=limits, grad = grad_fd)
m_fd.strategy = 0;
@btime migrad(m)
790.717 μs (706 allocations: 579.37 KiB)
Migrad | ||||
---|---|---|---|---|
FCN = -3.871e+04 | Nfcn = 115474 | |||
EDM = 2.01e-12 (Goal: 0.0002) | time = 0.1 sec | |||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | x0 | 10.0e3 | 0.1e3 | 0 | ||||
1 | x1 | -0.011 | 0.010 | |||||
2 | x2 | 0.999 | 0.007 | 0 |
x0 | x1 | x2 | |
---|---|---|---|
x0 | 9.91e+03 | 0.00143 | 0.00294 (0.004) |
x1 | 0.00143 | 0.0001 | -3.68e-07 (-0.005) |
x2 | 0.00294 (0.004) | -3.68e-07 (-0.005) | 5.1e-05 |
@btime migrad(m_fd)
721.547 μs (637 allocations: 521.48 KiB)
Migrad | ||||
---|---|---|---|---|
FCN = -3.871e+04 | Nfcn = 114537, Ngrad = 7 | |||
EDM = 1.5e-10 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | x0 | 10.0e3 | 0.1e3 | 0 | ||||
1 | x1 | -0.011 | 0.010 | |||||
2 | x2 | 0.999 | 0.008 | 0 |
x0 | x1 | x2 | |
---|---|---|---|
x0 | 9.94e+03 | 0.00354 | 0.094 (0.119) |
x1 | 0.00354 | 9.98e-05 | 2.17e-05 (0.275) |
x2 | 0.094 (0.119) | 2.17e-05 (0.275) | 6.24e-05 |