gistfile1.py

import numpy as np
import iminuit

def model(x, y, x_0, y_0, sigma_0, r_0, sigma_background):
    r = np.sqrt((x - x_0) ** 2 + (y - y_0) ** 2)
    return sigma_background + sigma_0 / (1 + r / r_0) ** 2

# Example "true" model parameters
x_0, y_0 = 0, 0 # Source position
r_0 = 10 # Source radius in pixels
sigma_0 = 100 # Source surface brightness at the center
sigma_background = 10 # Background surface brightness (flat distribution)
true_pars = (x_0, y_0, sigma_0, r_0, sigma_background)

# Simulate binned data on a grid with Poisson fluctuations
x_min, x_max = -100, 100
y_min, y_max = -80, 80
binsize = 0.5
xbin, ybin = np.meshgrid(np.arange(x_min, x_max, binsize),
                         np.arange(y_min, y_max, binsize))
counts = np.random.poisson(model(xbin, ybin, *true_pars))

def binned_nll(x_0, y_0, sigma_0, r_0, sigma_background):
    """Binned negative log likelihood fit statistic.
    Note that this is a function of the model parameters only,
    although to evaluate it the data is used."""
    mu = model(xbin, ybin, x_0, y_0, sigma_0, r_0, sigma_background)
    x = counts
    # This is the relevant part of -log(L) for Poisson likelihood L
    nll = np.where(mu > 0, -x * log(mu) + mu, 1)
    return nll.sum()

start_values = dict(x_0=x_0, y_0=y_0, sigma_0=sigma_0, r_0=r_0, sigma_background=sigma_background)
start_errors = dict(error_x_0=0.01, error_y_0=0.01, error_sigma_background=0.01)
fixed_parameters = dict(fix_r_0=True, fix_sigma_background=True)
pars = dict(start_values.items() + start_errors.items() + fixed_parameters.items())
minuit = iminuit.Minuit(binned_nll, pedantic=False, print_level=10, **pars)

print('*** Before MIGRAD *** ')
minuit.print_initial_param()
minuit.print_param()

print('*** Running MIGRAD *** ')
migrad = minuit.migrad(ncall=300)

print('*** After MIGRAD *** ')
minuit.print_param()

#minuit.fitarg