Index: trunk/npstat/swig/npstat_fitting.py
===================================================================
--- trunk/npstat/swig/npstat_fitting.py (revision 0)
+++ trunk/npstat/swig/npstat_fitting.py (revision 851)
@@ -0,0 +1,28 @@
+"""
+This module provides a number of useful wrapper classes and functions
+related to fitting probability densities with Minuit2. It requires
+"minuit2py" package.
+"""
+
+__author__ = "Igor Volobouev (i.volobouev@ttu.edu)"
+__version__ = "1.0"
+__date__ = "Oct 30 2022"
+
+import numpy as np
+import npstat as ns
+
+# The function to minimize should be represented by
+# a python callable. If the gradient is not calculated
+# explicitly, the callable should just return the value
+# of the function. Subsequently, the callable should
+# be used as an argument of the "Minuit2PySimpleFCN"
+# class constructor.
+class NegativeLogLikelihood:
+ def __init__(self, sample, distroConstructor, minDensity=2.2251e-308):
+ self.sample = sample
+ self.distroConstructor = distroConstructor
+ self.minDensity = minDensity
+ def __call__(self, *args):
+ distro = self.distroConstructor(*args)
+ densities = ns.scanDensity1D(distro, self.sample)
+ return -np.sum(np.log(np.maximum(densities, self.minDensity)))
Index: trunk/npstat/swig/npstat_plotting.py
===================================================================
--- trunk/npstat/swig/npstat_plotting.py (revision 850)
+++ trunk/npstat/swig/npstat_plotting.py (revision 851)
@@ -1,134 +1,144 @@
+"""
+This module provides a number of useful wrapper classes and functions
+related to plotting histograms and densities implemented in the "NPStat"
+nonparametric statistics package (Python 3 version).
+"""
+
+__author__ = "Igor Volobouev (i.volobouev@ttu.edu)"
+__version__ = "1.0"
+__date__ = "Oct 1 2022"
+
import npstat
import numpy as np
from npstat_utils import BoxND, arrayValueTypeAsNumpy, arrayNDFromNumpy
def _remove_keys(keys, dictionary):
for key in keys:
if key in dictionary:
del dictionary[key]
def colormeshHisto2D(ax, h, **options):
# We will flip the y axis in order to use the "normal" plotting
# convention instead of the matrix convention
assert h.dim() == 2, "Not a 2-d histogram"
xedges = h.axis(0).binEdges()
yedges = h.axis(1).binEdges(True)
bins = npstat.arrayNDToNumpy(h.binContents(), False)
xv, yv = np.meshgrid(xedges, yedges, indexing='ij')
# Get our options before passing **options to the pcolormesh
xlabel = options.get('xlabel', h.axis(0).label())
ylabel = options.get('ylabel', h.axis(1).label())
title = options.get('title', h.title())
_remove_keys(('xlabel', 'ylabel', 'title'), options)
mesh = ax.pcolormesh(xv, yv, np.flip(bins, 1), **options)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
return mesh
def colormeshDensity2D(ax, distro, bBox, shape, **options):
# We will flip the y axis in order to use the "normal" plotting
# convention instead of the matrix convention
assert distro.dim() == 2, "Not a 2-d distribution"
assert len(bBox) == 2, "Not a 2-d bounding box"
assert len(shape) == 2, "Not a 2-d shape"
xedges = np.linspace(bBox[0].min(), bBox[0].max(), shape[0]+1)
yedges = np.linspace(bBox[1].max(), bBox[1].min(), shape[1]+1)
# The following makes a numpy array
scan = npstat.scanDensityND(distro, bBox, shape)
xv, yv = np.meshgrid(xedges, yedges, indexing='ij')
xlabel = options.get('xlabel', 'X')
ylabel = options.get('ylabel', 'Y')
title = options.get('title', None)
_remove_keys(('xlabel', 'ylabel', 'title'), options)
mesh = ax.pcolormesh(xv, yv, np.flip(scan, 1), **options)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
if title is not None:
ax.set_title(title)
return mesh
def surfacePlotDensity2D(ax, distro, bBox, shape, **options):
# We will flip the y axis in order to use the "normal" plotting
# convention instead of the matrix convention
assert distro.dim() == 2, "Not a 2-d distribution"
assert len(bBox) == 2, "Not a 2-d bounding box"
assert len(shape) == 2, "Not a 2-d shape"
xedges = np.linspace(bBox[0].min(), bBox[0].max(), shape[0])
yedges = np.linspace(bBox[1].max(), bBox[1].min(), shape[1])
# The following makes a numpy array
scan = npstat.scanDensityND(distro, bBox, shape, True)
xv, yv = np.meshgrid(xedges, yedges, indexing='ij')
xlabel = options.get('xlabel', 'X')
ylabel = options.get('ylabel', 'Y')
zlabel = options.get('zlabel', 'Density')
title = options.get('title', None)
_remove_keys(('xlabel', 'ylabel', 'zlabel', 'title'), options)
surf = ax.plot_surface(xv, yv, np.flip(scan, 1), **options)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_zlabel(zlabel)
if title is not None:
ax.set_title(title)
return surf
def colormeshNumpyArray(ax, arr, bBox, **options):
assert arr.ndim == 2, "Not a 2-d array"
assert len(bBox) == 2, "Not a 2-d bounding box"
shape = arr.shape
xedges = np.linspace(bBox[0].min(), bBox[0].max(), shape[0]+1)
flip = options.get('flipY', False) or options.get('flipy', False)
if flip:
yedges = np.linspace(bBox[1].max(), bBox[1].min(), shape[1]+1)
else:
yedges = np.linspace(bBox[1].min(), bBox[1].max(), shape[1]+1)
xv, yv = np.meshgrid(xedges, yedges, indexing='ij')
xlabel = options.get('xlabel', 'X')
ylabel = options.get('ylabel', 'Y')
title = options.get('title', None)
_remove_keys(('xlabel', 'ylabel', 'title', 'flipY', 'flipy'), options)
if flip:
mesh = ax.pcolormesh(xv, yv, np.flip(arr, 1), **options)
else:
mesh = ax.pcolormesh(xv, yv, arr, **options)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
if title is not None:
ax.set_title(title)
return mesh
def plotHistoOutline1D(ax, h, **options):
xlabel = options.get('xlabel', h.axis(0).label())
ylabel = options.get('ylabel', h.accumulatedDataLabel())
title = options.get('title', h.title())
color = options.get('color', 'k')
width = h.axis(0).length()
xmin = options.get('xmin', h.axis(0).min() - 0.05*width)
xmax = options.get('xmax', h.axis(0).max() + 0.05*width)
x, y = npstat.histoOutline1D(h)
lowerbound = min(y)
upperbound = max(y)
if lowerbound > 0.0:
lowerbound = 0.0
yrange = upperbound - lowerbound
if yrange > 0.0:
up = upperbound + yrange*0.1
down = lowerbound - yrange*0.05
else:
up = 1.0
down = -0.25
ymin = options.get('ymin', down)
ymax = options.get('ymax', up)
_remove_keys(('xlabel', 'ylabel', 'title', 'color',
'xmin', 'xmax', 'ymin', 'ymax'), options)
ax.plot(x, y, color, **options)
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
def uncertaintyBandPolygon(x, ymin, ymax):
newx = np.concatenate((x, np.flip(x)))
newy = np.concatenate((ymin, np.flip(ymax)))
return newx, newy