{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Working in more than one dimension\n",
"Although in general the case of having a fit to 1D data is the most frequent, sometimes it is needed to perform fits in many dimensions. This is completely supported by MinKit, for both the unbinned and binned cases, although with some limitations. These limitations are imposed by the multidimensional nature of the problem, since calculating numerical integrals might lead to enormous arrays for the grids to use. As always, we start with a few imports, and we will define a function to help us to plot the results."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import minkit\n",
"import numpy as np\n",
"\n",
"def create_figure():\n",
"\n",
" fig = plt.figure(figsize=(10, 10))\n",
"\n",
" left, width = 0.12, 0.55\n",
" bottom, height = 0.12, 0.55\n",
" bottom_h = left_h = left + width + 0.02\n",
"\n",
" ax0 = plt.axes([left, bottom, width, height])\n",
" ax1 = plt.axes([left, bottom_h, width, 0.25], sharex=ax0)\n",
" ax2 = plt.axes([left_h, bottom, 0.25, height], sharey=ax0)\n",
"\n",
" plt.setp(ax1.get_xticklabels(), visible=False)\n",
" plt.setp(ax2.get_yticklabels(), visible=False)\n",
" \n",
" return fig, ax0, ax1, ax2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Unbinned case\n",
"Let's create an unbinned data set in two dimensions from two different Gaussian distributions with shared mean and different standard deviation, do a fit to it, and plot the results."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" FCN = 694289.713121555 \n",
" TOTAL NCALL = 81 \n",
" NCALLS = 81 \n",
" \n",
" \n",
" EDM = 9.269661053992134e-05 \n",
" GOAL EDM = 1e-05 \n",
" \n",
" UP = 1.0 \n",
" \n",
"
\n",
"\n",
" \n",
" Valid \n",
" Valid Param \n",
" Accurate Covar \n",
" PosDef \n",
" Made PosDef \n",
" \n",
" \n",
" True \n",
" True \n",
" True \n",
" True \n",
" False \n",
" \n",
" \n",
" Hesse Fail \n",
" HasCov \n",
" Above EDM \n",
" Reach calllim \n",
" \n",
" \n",
" False \n",
" True \n",
" False \n",
" False \n",
" \n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
" \n",
" + Name \n",
" Value \n",
" Hesse Error \n",
" Minos Error- \n",
" Minos Error+ \n",
" Limit- \n",
" Limit+ \n",
" Fixed? \n",
" \n",
" \n",
" 0 \n",
" c \n",
" 0.00187518 \n",
" 0.00285409 \n",
" -5 \n",
" 5 \n",
" No \n",
" \n",
" \n",
" 1 \n",
" sx \n",
" 1.99467 \n",
" 0.00544549 \n",
" 1 \n",
" 3 \n",
" No \n",
" \n",
" \n",
" 2 \n",
" sy \n",
" 0.999087 \n",
" 0.00223423 \n",
" 0.5 \n",
" 3 \n",
" No \n",
" \n",
"
\n",
"\n",
"\n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Define the model\n",
"x = minkit.Parameter('x', bounds=(-5, +5))\n",
"y = minkit.Parameter('y', bounds=(-5, +5))\n",
"c = minkit.Parameter('c', 0, bounds=(-5, +5))\n",
"\n",
"sx = minkit.Parameter('sx', 2, bounds=(1, 3))\n",
"sy = minkit.Parameter('sy', 1, bounds=(0.5, 3))\n",
"\n",
"gx = minkit.Gaussian('gx', x, c, sx)\n",
"gy = minkit.Gaussian('gy', y, c, sy)\n",
"\n",
"pdf = minkit.ProdPDFs('pdf', [gx, gy])\n",
"\n",
"# Generate data and fit it\n",
"data = pdf.generate(100000)\n",
"\n",
"with minkit.minimizer('uml', pdf, data) as minimizer:\n",
" minimizer.migrad()\n",
"\n",
"# Plot the results\n",
"values, (xedges, yedges) = minkit.data_plotting_arrays(data, bins=100)\n",
"\n",
"cx = 0.5 * (xedges[1:] + xedges[:-1])\n",
"cy = 0.5 * (yedges[1:] + yedges[:-1])\n",
"\n",
"vx, vy = tuple(a.flatten() for a in np.meshgrid(cx, cy))\n",
"\n",
"fig, ax0, ax1, ax2 = create_figure()\n",
"\n",
"# Plot the data\n",
"ax0.hist2d(vx, vy, bins=(xedges, yedges), weights=values);\n",
"\n",
"# Calculate the projections of the data and plot them\n",
"xproj, _ = minkit.data_plotting_arrays(data, bins=100, projection='x')\n",
"yproj, _ = minkit.data_plotting_arrays(data, bins=100, projection='y')\n",
"\n",
"ax1.hist(cx, bins=xedges, weights=xproj, histtype='step', color='k');\n",
"ax2.hist(cy, bins=yedges, weights=yproj, histtype='step', orientation='horizontal', color='k');\n",
"\n",
"# Calculate the values of the PDF and plot them\n",
"gxc, xpdf = minkit.pdf_plotting_arrays(pdf, values, (xedges, yedges), projection='x')\n",
"gyc, ypdf = minkit.pdf_plotting_arrays(pdf, values, (xedges, yedges), projection='y')\n",
"\n",
"ax1.plot(gxc, xpdf, color='b');\n",
"ax2.plot(ypdf, gyc, color='b');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The amount of code needed to generate and fit is very similar to that of the 1D case. In order to plot the results, it is needed to work a bit more if one wants to get the projections of the PDF in each dimension, since it comes necessary to sum the contribution in the other. By using *numpy.histogramdd* one can extend the previous code to the n-dimensional case, although computation might turn slower due to the curse of dimensionality."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Binned case\n",
"Similarly to what has been done in the previous section, one can do the same thing on binned data sets."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" FCN = 7115.751392859857 \n",
" TOTAL NCALL = 34 \n",
" NCALLS = 34 \n",
" \n",
" \n",
" EDM = 2.6395486866473276e-09 \n",
" GOAL EDM = 1e-05 \n",
" \n",
" UP = 1.0 \n",
" \n",
"
\n",
"\n",
" \n",
" Valid \n",
" Valid Param \n",
" Accurate Covar \n",
" PosDef \n",
" Made PosDef \n",
" \n",
" \n",
" True \n",
" True \n",
" True \n",
" True \n",
" False \n",
" \n",
" \n",
" Hesse Fail \n",
" HasCov \n",
" Above EDM \n",
" Reach calllim \n",
" \n",
" \n",
" False \n",
" True \n",
" False \n",
" False \n",
" \n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
" \n",
" + Name \n",
" Value \n",
" Hesse Error \n",
" Minos Error- \n",
" Minos Error+ \n",
" Limit- \n",
" Limit+ \n",
" Fixed? \n",
" \n",
" \n",
" 0 \n",
" c \n",
" 0.00191893 \n",
" 0.00285086 \n",
" -5 \n",
" 5 \n",
" No \n",
" \n",
" \n",
" 1 \n",
" sx \n",
" 1.99482 \n",
" 0.00542073 \n",
" 1 \n",
" 3 \n",
" No \n",
" \n",
" \n",
" 2 \n",
" sy \n",
" 0.999001 \n",
" 0.00223591 \n",
" 0.5 \n",
" 3 \n",
" No \n",
" \n",
"
\n",
"\n",
"\n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Transform data to binned, and do a fit to it\n",
"binned_data = data.make_binned(bins=(100, 100))\n",
"\n",
"xedges, yedges = binned_data['x'].as_ndarray(), binned_data['y'].as_ndarray()\n",
"\n",
"values = binned_data.values.as_ndarray()\n",
"\n",
"with minkit.minimizer('bml', pdf, binned_data) as minimizer:\n",
" minimizer.migrad()\n",
"\n",
"# Plot the results\n",
"fig, ax0, ax1, ax2 = create_figure()\n",
"\n",
"data_values, (xedges, yedges) = minkit.data_plotting_arrays(binned_data)\n",
"\n",
"cx = 0.5 * (xedges[1:] + xedges[:-1])\n",
"cy = 0.5 * (yedges[1:] + yedges[:-1])\n",
"\n",
"vx, vy = tuple(a.flatten() for a in np.meshgrid(cx, cy))\n",
"\n",
"ax0.hist2d(vx, vy, bins=(xedges, yedges), weights=data_values);\n",
"\n",
"# Calculate the projections of the data and plot them\n",
"xproj, _ = minkit.data_plotting_arrays(binned_data, projection='x')\n",
"yproj, _ = minkit.data_plotting_arrays(binned_data, projection='y')\n",
"\n",
"cx = 0.5 * (xedges[1:] + xedges[:-1])\n",
"cy = 0.5 * (yedges[1:] + yedges[:-1])\n",
"\n",
"ax1.hist(cx, bins=xedges, weights=xproj, histtype='step', color='k');\n",
"ax2.hist(cy, bins=yedges, weights=yproj, histtype='step', orientation='horizontal', color='k');\n",
"\n",
"# Calculate the values of the PDF and plot them\n",
"gxc, xpdf = minkit.pdf_plotting_arrays(pdf, values, (xedges, yedges), projection='x')\n",
"gyc, ypdf = minkit.pdf_plotting_arrays(pdf, values, (xedges, yedges), projection='y')\n",
"\n",
"ax1.plot(gxc, xpdf, color='b');\n",
"ax2.plot(ypdf, gyc, color='b');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As it can be clearly seen, the result is similar to that of the unbinned case."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}