diff --git a/docs/examples/notebooks/custom_modifiers.ipynb b/docs/examples/notebooks/custom_modifiers.ipynb index 6a91e4f5ca..118b9d5a94 100644 --- a/docs/examples/notebooks/custom_modifiers.ipynb +++ b/docs/examples/notebooks/custom_modifiers.ipynb @@ -1 +1,184 @@ -{"cells":[{"cell_type":"code","execution_count":1,"metadata":{},"outputs":[],"source":["import pyhf\n","import matplotlib.pyplot as plt\n","import numpy as np\n","import scipy.stats"]},{"cell_type":"code","execution_count":21,"metadata":{},"outputs":[],"source":["class custom_modifier(object):\n"," op_code = 'addition'\n"," def __init__(self,name = 'name'):\n"," self.required_parsets = {\n"," 'mean': {\n"," 'paramset_type': pyhf.parameters.paramsets.unconstrained,\n"," 'n_parameters': 1,\n"," 'is_constrained': False,\n"," 'is_shared': True,\n"," 'inits': (0.0,),\n"," 'bounds': ((-5, 5),),\n"," 'fixed': False\n"," }}\n"," def apply(self,pars):\n"," base = np.zeros((1,len(m.config.samples),1,sum(m.config.channel_nbins.values())))\n","\n"," bins = np.linspace(-5,5,20+1)\n"," mean = pars[self.config.par_slice('mean')][0]\n"," yields = 100*(scipy.stats.norm(loc = mean).cdf(bins[1:]) - scipy.stats.norm(loc = mean).cdf(bins[:-1]))\n"," base[0,self.config.samples.index('signal'),0,:] = yields\n"," return base\n","\n","\n","\n","m = pyhf.Model(\n","{'channels': [{'name': 'singlechannel',\n"," 'samples': [{'name': 'signal',\n"," 'data': [0]*20,\n"," 'modifiers': [{'name': 'mu', 'type': 'normfactor', 'data': None}]},\n"," {'name': 'background',\n"," 'data': [300]*20,\n"," 'modifiers': []}]}]},\n"," custom_modifiers=[custom_modifier()]\n",")\n","bp = pyhf.tensorlib.astensor(m.config.suggested_init())\n","bp[m.config.poi_index] = 5\n","bp[m.config.par_slice('mean')] = [3.0]\n","d = m.make_pdf(bp).sample()"]},{"cell_type":"code","execution_count":28,"metadata":{},"outputs":[{"output_type":"stream","name":"stdout","text":["[3.06386254]\n[5.17062475]\n"]}],"source":["bestfit = pyhf.infer.mle.fit(d,m)\n","print(bestfit[m.config.par_slice('mean')])\n","print(bestfit[m.config.par_slice('mu')])\n"]},{"cell_type":"code","execution_count":30,"metadata":{},"outputs":[{"output_type":"execute_result","data":{"text/plain":[""]},"metadata":{},"execution_count":30},{"output_type":"display_data","data":{"text/plain":"
","image/svg+xml":"\n\n\n\n \n \n \n \n 2020-11-20T14:04:05.443973\n image/svg+xml\n \n \n Matplotlib v3.3.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n","image/png":"iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAb4ElEQVR4nO3df7BU5Z3n8feHC4HRkIjCKnIRiDETiSmQXH8kGmF1J4JlgSYZg2slmHFzR9CqpGYSo2PVxJ0sVf5K3HVrIIWrK44k6piIxCIhapQkldUIDhLAOKKBcCmQHxolGh3B7/7Rz8X20n1v//5xz+dV1XXPec5z+nz73O5vn37Oc56jiMDMzAa/Ic0OwMzMGsMJ38wsI5zwzcwywgnfzCwjnPDNzDJiaLMDABg9enRMnDix2WGYmbWVtWvX7omIMaXWb4mEP3HiRNasWdPsMMzM2oqkreXUd5OOmVlGOOGbmWWEE76ZWUa0RBt+IW+//TY9PT28+eabzQ6lrkaMGEFnZyfDhg1rdihmNsi1bMLv6elh5MiRTJw4EUnNDqcuIoK9e/fS09PDpEmTmh2OmQ1yLduk8+abb3LUUUcN2mQPIImjjjpq0P+KMbPW0LIJHxjUyb5XFl6jmbWGlk74ZmZWO074Zbjuuuu4+eabiy5fvnw5mzZtamBEZmalc8KvISd8s/qaMWMGM2bMaHYYbWvQJPxly5YxceJEhgwZwsSJE1m2bFlNnnfhwoV85CMf4cwzz+S5554D4LbbbuOUU05hypQpfO5zn+ONN97g17/+NStWrOAb3/gGU6dO5YUXXihYz8ysWQZFwl+2bBnd3d1s3bqViGDr1q10d3dXnfTXrl3LPffcw7p161i5ciVPPfUUAJ/97Gd56qmneOaZZzjxxBO5/fbb+dSnPsXs2bO56aabWLduHccff3zBemZmzTIoEv611157yNHzG2+8wbXXXlvV8/7yl7/kwgsv5LDDDuMDH/gAs2fPBmDDhg18+tOf5uMf/zjLli1j48aNBdcvtZ6ZWSO07IVX5fjDH/5QVnm1Lr30UpYvX86UKVO48847efzxx6uqZ2bWCCUf4UvqkPRvkh5K85MkPSlps6R7Jb0vlQ9P85vT8ol1iv2g4447rqzyUp111lksX76cP//5z+zbt48f//jHAOzbt4+xY8fy9ttvv6fZaOTIkezbt+/gfLF6Zlnlk67NVU6TzleBZ/PmbwBuiYgPA68Al6Xyy4BXUvktqV5dLVy4kMMOO+w9ZYcddhgLFy6s6nmnTZvGF77wBaZMmcKsWbM45ZRTAPj2t7/NaaedxhlnnMFHP/rRg/Xnzp3LTTfdxMknn8wLL7xQtJ6ZWVNExIAPoBN4FDgbeAgQsAcYmpZ/EliVplcBn0zTQ1M99ff8n/jEJ6KvTZs2HVLWn7vvvjsmTJgQkmLChAlx9913l7V+M5X7Ws3a1fTp02P69OkVrTt//vwAAoiOjo6YP39+bYNrQ8CaKCGH9z5KbcP/n8BVwMg0fxTwx4jYn+Z7gHFpehywLX2Z7Jf0aqq/p8zvorJccsklXHLJJfXchJk1yNjO49i5fVvR5QcOHGDx4sUsXrz4YNkx48azo6c+5+0GiwETvqTzgV0RsVbSjFptWFI30A3Vt7Wb2eCyc/s2JnzzoYPzW2+cDfHOoRU1hAlXrcjVueH8RoXXtko5wj8DmC3pPGAE8AHgfwFHSBqajvI7ge2p/nZgPNAjaSjwQWBv3yeNiCXAEoCurq6o9oWYWWvp7yi90KCB/R6hF0r2/ZVbQQMm/Ii4BrgGIB3hfz0iLpH0r8DngXuAecCDaZUVaf7/peU/T21NZpYhfY/SAXZ+/2oAjvmv1x9Sv98jdA0peoRvpatmb30T+DtJm8m10fdeRno7cFQq/zvg6upCNLNaalbXyD2rFvHWtg28tW0DW2+czZ5Vi0pe9/ApM8sqt8LKuvAqIh4HHk/TLwKnFqjzJvDXNYjNzAaJPasW8fq6le8WxDsH50efu2DA9XvrHHwODeHwKTNLWtfe1Ta/h8Z2Hoekmj3Gdg58ovjWW2/lxBNPZNSoUVx/fe4nqEfENCvf68/8tKzyQkafu4Dh409i+PiTmHDVioYm+8FywVjbDK1QqD2wGqWc0V+0aBGPPPIInZ2dB8uWL1/O+eefz+TJk2sWi1k76U18ZQ0V4pOuLaFtjvAb7fLLL+fFF19k1qxZ3HLLLVx55ZUFh0A2sxIUO7nawJOug+UovRpO+EV873vf49hjj+Wxxx5j1KhRAAWHQDazgfmka2tomyYdM2tfjT7pWtNrAAYRJ3yzDFmwYAGrV68GYOjQoXR3d7NoUendI6sx+twF7N+bS6qF+uGXotT1anoNwCDiJp0y9R0C2axdLFiw4D1jz/SOR7Nggbs2ZkXbHOEfM258Tb+Fjxk3vqL15s6dy1e+8hVuvfVW7r//frfjW9tYsmRJ0fJSj/Kb+QuhWQbTa26bhN+M9rUtW7YAuTtXXXrppQCcccYZ7odvbenAgQNllUP/beHtNGJl71W+kBuIrdTzB8V+FQFtmfTdpGOWER0dHWWVw7tt4RO++VC/XSt76/Q3pHGzFLvKt5ShHfr7VdSOnPDNMqK7u7us8kO06cVT1VzlW8mvolbW0k06EVGwC9Vg4oFErVF6myB6myQ6OjrKa4+uwYiVlfbOqUoVX1QdHR0Fk3t/v4paWcse4Y8YMYK9e/cO6oQYEezdu5cRI0Y0OxTLiEWLFjF9+nSmT5/O/v37y2qHbtuLp6q4yrfqX0UtpmWP8Ds7O+np6WH37t3NDqWuRowY8Z6xesxaVbuOWHn4lJnvbcPPKx9I1b+KWkzLJvxhw4YxadKkZodhZnlqcfFUo1X7RbVo0aKDPfPKGjAuT0UDztVByyZ8M7NaaccvqnoYsBFL0ghJv5H0jKSNkv57Kr9T0u8lrUuPqalckm6VtFnSeknT6vwazMysBKUc4b8FnB0Rf5I0DPiVpJ+kZd+IiPv71J8FnJAepwGL018zM2uiUm5iHsCf0uyw9Oiv68wc4K603hOSjpA0NiJ2VB2tmVWt2e3I1jwldcuU1CFpHbALeDginkyLFqZmm1skDU9l44D8y+16Ulnf5+yWtEbSmsHeE8fM2keh26muXr2a1atXV3y71FZR0knbiDgATJV0BPCApJOAa4CdwPuAJcA3gX8qdcMRsSStR1dX1+DtbG9WRKv03CjXYD/pWe7tVAca1LGVBl8rq5dORPxR0mPAzIi4ORW/Jen/Al9P89uB/KEoO1OZmVnTNOOLqtUGXyull86YdGSPpL8A/gr4naSxqUzABcCGtMoK4Eupt87pwKtuv2+MrN6zM6uvuxSFmicGerRTE0Wra7XB10o5wh8LLJXUQe4L4r6IeEjSzyWNAQSsAy5P9VcC5wGbgTeAL9c8ajMrSbnNE5Cduz81QqsNvlZKL531wMkFys8uUj+AK6oPzcysvbXa4GstO3iamVm7a7XB15zwzZqgt+fG6tWrGTp0qO8rO0gtWrSI+fPnH5zv6Ohg/vz57dFLx8yq12o9N6y+ajH4Wq34CL+Aduz1kdUjxnZ83a3Wc8Oywwl/ECh2xNgOya8a7fq6W63nhmWHE/4gkNUjxnZ93ZXcTNysFpzwB4GsHjE2+3VX2vTXaj03LDt80nYQaLW+vo3Srq+7FrfNa9dxeAa7sZ3HsXP7tqLLcwMTvNcx48azo+cP9QzroEGb8LP0geju7n5PW3Z+eaM0Y3+3wuuuVCv13LDaafUrmwdtwq9Us0e2qyRxDrYbLZcqK6+7v6PGZh8xWntxws/Tzv2js3rEmIXXXeiocef3rwYKjwDpsXCsGJ+0zdOuvT7MzErhhJ+nFr0+2vGiLWhu3O247Xa84MtsUDbpVNoO3669PqyxatX0V2kT1J5Vi3hrW+72E1tvnM3hU2Yy+lx/4djABt0RfjVXX7p/dGWydrTbzKa/PasW8fq6le8WxDu8vm4le1a19jkmaw2DLuFX82Fs9sh21SbOxx9/vOEnLltheINGv+5Kmv7KvfNUsbtOvf7MT8sqN8s3YJOOpBHAL4Dhqf79EfEtSZOAe4CjgLXAFyPiPyQNB+4CPgHsBb4QEVvqFP8hqm2Hb1avj3btIdTfF2wrx12NSpr+anZj7HinvHKzPKUc4b8FnB0RU4CpwMx0r9obgFsi4sPAK8Blqf5lwCup/JZUr2GaOU5JNUfozWwmqCbuZg9v0AxNbfpTkY9ssXKzPAO+SyLnT2l2WHoEcDZwfypfSu5G5gBz0jxp+TkqdHVInTTrw1ht00azEme1cVf7BdvM9v9ytp3fJFPo6l7IXQBW7xuBHz5lZlnlZvlK6qWTbmC+Fvgw8M/AC8AfI2J/qtIDjEvT44BtABGxX9Kr5Jp99vR5zm6gG+C44yr/cAw0dkWvxYsXH/yg5l+JWM1VjKVsu9h2y4m9N45yt12ruCtZ/8CBA2XH3ftlU812+6rVtlvl4qfe3jgHT9xqiHvpWMlKSvgRcQCYKukI4AHgo9VuOCKWAEsAurq6otLnKdY2WuqHsZoPcv66/X3Ai9XJX/+Q3hfJ4VPPO/hhHiju/lQT90Dbfk/sBRJQ0W3fOLtw27OGMOGqFQNutxS12nYrGX3uAvbvzX0RFXqPmhVTVj/8iPijpMeATwJHSBqajvI7ge2p2nZgPNAjaSjwQXInb9tGRR8iDSmaQErRtCO3KuOGKhJQM09A+uSnZdCAn2pJY9KRPZL+Avgr4FngMeDzqdo84ME0vSLNk5b/PCIqPoJvF7VoWx197gKGjz+J4eNPYsJVKxryM72pbcLNPAHpk5+WQaUc4Y8FlqZ2/CHAfRHxkKRNwD2S/gfwb8Dtqf7twL9I2gy8DMytQ9wtp13bVpsZ9+FTZhZuxmrAl00ttu3mFGs3Ayb8iFgPnFyg/EXg1ALlbwJ/XZPo2ky7tq02K+5mftm06xd0r3Z6f1nrGJRj6YA/EI1W6f5u5pdku35Bm1Vq0Cb8duXEY2b14jNUZmYZ4YRvZpYRbtKpsXZtkmlm3Fndtlmj+QjfzCwjnPDNzDLCCd/MLCOc8M3MMsIJ38wsI5zwzcwywgnfzCwjnPDNzDLCCd/MLCOc8M3MMsIJ38wsI0q5xeF4SY9J2iRpo6SvpvLrJG2XtC49zstb5xpJmyU9J+ncer4AMzMrTSmDp+0H/j4inpY0Elgr6eG07JaIuDm/sqTJ5G5r+DHgWOARSR+JiAO1DNzMzMoz4BF+ROyIiKfT9D5yNzAf188qc4B7IuKtiPg9sJkCt0I0M7PGKqsNX9JEcve3fTIVXSlpvaQ7JI1KZeOAbXmr9VDgC0JSt6Q1ktbs3r27/MjNzKwsJSd8Se8Hfgh8LSJeAxYDxwNTgR3Ad8rZcEQsiYiuiOgaM2ZMOauamVkFSkr4koaRS/bLIuJHABHxUkQciIh3gNt4t9lmOzA+b/XOVGZmZk1USi8dAbcDz0bEd/PKx+ZVuxDYkKZXAHMlDZc0CTgB+E3tQjYzs0qU0kvnDOCLwG8lrUtl/wBcLGkqEMAW4G8BImKjpPuATeR6+FzhHjpmZs03YMKPiF8BKrBoZT/rLAQWVhGXmZnVmK+0NTPLCCd8M7OMcMI3M8sIJ3wzs4xwwjczywgnfDOzjHDCNzPLCCd8M7OMcMI3M8sIJ3wzs4xwwjczywgnfDOzjHDCNzPLCCd8M7OMcMI3M8sIJ3wzs4wo5RaH4yU9JmmTpI2SvprKj5T0sKTn099RqVySbpW0WdJ6SdPq/SLMzGxgpRzh7wf+PiImA6cDV0iaDFwNPBoRJwCPpnmAWeTuY3sC0A0srnnUZmZWtgETfkTsiIin0/Q+4FlgHDAHWJqqLQUuSNNzgLsi5wngiD43PDczsyYoqw1f0kTgZOBJ4OiI2JEW7QSOTtPjgG15q/Wksr7P1S1pjaQ1u3fvLjduMzMrU8kJX9L7gR8CX4uI1/KXRUQAUc6GI2JJRHRFRNeYMWPKWdXMzCpQUsKXNIxcsl8WET9KxS/1NtWkv7tS+XZgfN7qnanMzMyaqJReOgJuB56NiO/mLVoBzEvT84AH88q/lHrrnA68mtf0Y2ZmTTK0hDpnAF8EfitpXSr7B+B64D5JlwFbgYvSspXAecBm4A3gy7UM2MzMKjNgwo+IXwEqsvicAvUDuKLKuMzMrMZ8pa2ZWUY44ZuZZYQTvplZRjjhm5llhBO+mVlGOOGbmWWEE76ZWUY44ZuZZYQTvplZRjjhm5llhBO+mVlGOOGbmWWEE76ZWUY44ZuZZYQTvplZRpRyx6s7JO2StCGv7DpJ2yWtS4/z8pZdI2mzpOcknVuvwM3MrDylHOHfCcwsUH5LRExNj5UAkiYDc4GPpXUWSeqoVbBmZla5ARN+RPwCeLnE55sD3BMRb0XE78nd5vDUKuIzM7MaqaYN/0pJ61OTz6hUNg7YllenJ5UdQlK3pDWS1uzevbuKMMzMrBSVJvzFwPHAVGAH8J1ynyAilkREV0R0jRkzpsIwzMysVBUl/Ih4KSIORMQ7wG2822yzHRifV7UzlZmZWZNVlPAljc2bvRDo7cGzApgrabikScAJwG+qC9HMzGph6EAVJP0AmAGMltQDfAuYIWkqEMAW4G8BImKjpPuATcB+4IqIOFCXyM3MrCwDJvyIuLhA8e391F8ILKwmKDMzqz1faWtmlhFO+GZmGeGEb2aWEU74ZmYZ4YRvZpYRTvhmZhnhhG9mlhFO+GZmGeGEb2aWEU74ZmYZ4YRvZpYRTvhmZhnhhG9mlhFO+GZmGeGEb2aWEU74ZmYZMWDCl3SHpF2SNuSVHSnpYUnPp7+jUrkk3Spps6T1kqbVM3gzMytdKUf4dwIz+5RdDTwaEScAj6Z5gFnk7mN7AtANLK5NmGZmVq0BE35E/AJ4uU/xHGBpml4KXJBXflfkPAEc0eeG52Zm1iSVtuEfHRE70vRO4Og0PQ7YllevJ5UdQlK3pDWS1uzevbvCMMzMrFRVn7SNiACigvWWRERXRHSNGTOm2jDMzGwAlSb8l3qbatLfXal8OzA+r15nKjMzsyarNOGvAOal6XnAg3nlX0q9dU4HXs1r+jEzsyYaOlAFST8AZgCjJfUA3wKuB+6TdBmwFbgoVV8JnAdsBt4AvlyHmM3MrAIDJvyIuLjIonMK1A3gimqDMjOz2vOVtmZmGeGEb2aWEU74ZmYZ4YRvZpYRTvhmZhnhhG9mlhFO+GZmGeGEb2aWEU74ZmYZ4YRvZpYRTvhmZhnhhG9mlhFO+GZmGeGEb2aWEU74ZmYZ4YRvZpYRA94ApT+StgD7gAPA/ojoknQkcC8wEdgCXBQRr1QXppmZVasWR/j/OSKmRkRXmr8aeDQiTgAeTfNmZtZk9WjSmQMsTdNLgQvqsA0zMytTtQk/gJ9JWiupO5UdHRE70vRO4OhCK0rqlrRG0prdu3dXGYaZmQ2kqjZ84MyI2C7pPwEPS/pd/sKICElRaMWIWAIsAejq6ipYx8zMaqeqI/yI2J7+7gIeAE4FXpI0FiD93VVtkGZmVr2KE76kwyWN7J0GPgNsAFYA81K1ecCD1QZpZmbVq6ZJ52jgAUm9z/P9iPippKeA+yRdBmwFLqo+TDMzq1bFCT8iXgSmFCjfC5xTTVBmZlZ7vtLWzCwjnPDNzDLCCd/MLCOc8M3MMsIJ38wsI5zwzcwywgnfzCwjnPDNzDLCCd/MLCOc8M3MMsIJ38wsI5zwzcwywgnfzCwjnPDNzDLCCd/MLCOc8M3MMqJuCV/STEnPSdos6ep6bcfMzEpTl4QvqQP4Z2AWMBm4WNLkemzLzMxKU68j/FOBzRHxYkT8B3APMKdO2zIzsxIoImr/pNLngZkR8d/S/BeB0yLiyrw63UB3mv1L4LmaBwKjgT11eN5qOa7ytWpsrRoXtG5srRoXtG5sxeKaEBFjSn2Sim9iXq2IWAIsqec2JK2JiK56bqMSjqt8rRpbq8YFrRtbq8YFrRtbreKqV5POdmB83nxnKjMzsyapV8J/CjhB0iRJ7wPmAivqtC0zMytBXZp0ImK/pCuBVUAHcEdEbKzHtgZQ1yajKjiu8rVqbK0aF7RubK0aF7RubDWJqy4nbc3MrPX4Slszs4xwwjczy4i2T/gDDeEgabike9PyJyVNbFBc4yU9JmmTpI2SvlqgzgxJr0palx7/2KDYtkj6bdrmmgLLJenWtM/WS5rWoLj+Mm9frJP0mqSv9anTkH0m6Q5JuyRtyCs7UtLDkp5Pf0cVWXdeqvO8pHkNiu0mSb9L/68HJB1RZN1+//d1iOs6Sdvz/l/nFVm3rkOxFInt3ry4tkhaV2Tdeu6zgnmibu+1iGjbB7kTwi8AHwLeBzwDTO5TZwHwvTQ9F7i3QbGNBaal6ZHAvxeIbQbwUBP22xZgdD/LzwN+Agg4HXiySf/bneQuLGn4PgPOAqYBG/LKbgSuTtNXAzcUWO9I4MX0d1SaHtWA2D4DDE3TNxSKrZT/fR3iug74egn/634/x/WIrc/y7wD/2IR9VjBP1Ou91u5H+KUM4TAHWJqm7wfOkaR6BxYROyLi6TS9D3gWGFfv7dbIHOCuyHkCOELS2AbHcA7wQkRsbfB2AYiIXwAv9ynOfy8tBS4osOq5wMMR8XJEvAI8DMysd2wR8bOI2J9mnyB37UtDFdlnpaj7UCz9xZbywUXAD2q5zVL0kyfq8l5r94Q/DtiWN9/DoUn1YJ30gXgVOKoh0SWpGelk4MkCiz8p6RlJP5H0sQaFFMDPJK1VboiLvkrZr/U2l+IfwGbsM4CjI2JHmt4JHF2gTivsu78h9wutkIH+9/VwZWpquqNI00Sz99mngZci4vkiyxuyz/rkibq819o94bc8Se8Hfgh8LSJe67P4aXJNFlOA/w0sb1BYZ0bENHKjmV4h6awGbbckyl2sNxv41wKLm7XP3iNyv6lbrk+zpGuB/cCyIlUa/b9fDBwPTAV2kGs6aTUX0//Rfd33WX95opbvtXZP+KUM4XCwjqShwAeBvY0ITtIwcv/EZRHxo77LI+K1iPhTml4JDJM0ut5xRcT29HcX8AC5n9T5mj00xizg6Yh4qe+CZu2z5KXepq30d1eBOk3bd5IuBc4HLklJ4hAl/O9rKiJeiogDEfEOcFuR7TVznw0FPgvcW6xOvfdZkTxRl/dauyf8UoZwWAH0nr3+PPDzYh+GWkrtgrcDz0bEd4vUOab3fIKkU8n9P+r6ZSTpcEkje6fJnezb0KfaCuBLyjkdeDXv52UjFD3iasY+y5P/XpoHPFigzirgM5JGpeaLz6SyupI0E7gKmB0RbxSpU8r/vtZx5Z/7ubDI9po5FMt/AX4XET2FFtZ7n/WTJ+rzXqvHmedGPsj1KPl3cmf5r01l/0TujQ8wglzTwGbgN8CHGhTXmeR+hq0H1qXHecDlwOWpzpXARnK9Ep4APtWAuD6UtvdM2nbvPsuPS+RuYPMC8Fugq4H/z8PJJfAP5pU1fJ+R+8LZAbxNrm30MnLnfh4FngceAY5MdbuA/5O37t+k99tm4MsNim0zufbc3vdab8+0Y4GV/f3v6xzXv6T30HpySWxs37jS/CGf43rHlsrv7H1v5dVt5D4rlifq8l7z0ApmZhnR7k06ZmZWIid8M7OMcMI3M8sIJ3wzs4xwwjczywgnfDOzjHDCNzPLiP8PhAA05t/oWvAAAAAASUVORK5CYII=\n"},"metadata":{"needs_background":"light"}}],"source":["plt.bar(np.arange(20),m.expected_actualdata(bestfit), alpha = 1.0, facecolor = None, edgecolor = 'k', label = 'fit')\n","plt.scatter(np.arange(20),d[:m.config.nmaindata], alpha = 1.0, marker = 'o', c='k', label = 'data', zorder=99)\n","plt.errorbar(np.arange(20),d[:m.config.nmaindata],yerr=np.sqrt(d[:m.config.nmaindata]), marker='o', c='k', linestyle = '')\n","plt.legend()\n"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["s"]}],"nbformat":4,"nbformat_minor":2,"metadata":{"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.7.2-final"},"orig_nbformat":2,"kernelspec":{"name":"python3","display_name":"Python 3"}}} \ No newline at end of file +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import pyhf\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import scipy.stats" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "class custom_modifier:\n", + " op_code = 'addition'\n", + "\n", + " def __init__(self, name='name'):\n", + " self.required_parsets = {\n", + " 'mean': {\n", + " 'paramset_type': pyhf.parameters.paramsets.unconstrained,\n", + " 'n_parameters': 1,\n", + " 'is_constrained': False,\n", + " 'is_shared': True,\n", + " 'inits': (0.0,),\n", + " 'bounds': ((-5, 5),),\n", + " 'fixed': False,\n", + " }\n", + " }\n", + "\n", + " def apply(self, pars):\n", + " base = np.zeros(\n", + " (1, len(m.config.samples), 1, sum(m.config.channel_nbins.values()))\n", + " )\n", + "\n", + " bins = np.linspace(-5, 5, 20 + 1)\n", + " mean = pars[self.config.par_slice('mean')][0]\n", + " yields = 100 * (\n", + " scipy.stats.norm(loc=mean).cdf(bins[1:])\n", + " - scipy.stats.norm(loc=mean).cdf(bins[:-1])\n", + " )\n", + " base[0, self.config.samples.index('signal'), 0, :] = yields\n", + " return base\n", + "\n", + "\n", + "m = pyhf.Model(\n", + " {\n", + " 'channels': [\n", + " {\n", + " 'name': 'singlechannel',\n", + " 'samples': [\n", + " {\n", + " 'name': 'signal',\n", + " 'data': [0] * 20,\n", + " 'modifiers': [\n", + " {'name': 'mu', 'type': 'normfactor', 'data': None}\n", + " ],\n", + " },\n", + " {'name': 'background', 'data': [300] * 20, 'modifiers': []},\n", + " ],\n", + " }\n", + " ]\n", + " },\n", + " custom_modifiers=[custom_modifier()],\n", + ")\n", + "bp = pyhf.tensorlib.astensor(m.config.suggested_init())\n", + "bp[m.config.poi_index] = 5\n", + "bp[m.config.par_slice('mean')] = [3.0]\n", + "d = m.make_pdf(bp).sample()" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "[3.06386254]\n[5.17062475]\n" + ] + } + ], + "source": [ + "bestfit = pyhf.infer.mle.fit(d, m)\n", + "print(bestfit[m.config.par_slice('mean')])\n", + "print(bestfit[m.config.par_slice('mu')])" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "" + ] + }, + "metadata": {}, + "execution_count": 30 + }, + { + "output_type": "display_data", + "data": { + "text/plain": "
", + "image/svg+xml": "\n\n\n\n \n \n \n \n 2020-11-20T14:04:05.443973\n image/svg+xml\n \n \n Matplotlib v3.3.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAb4ElEQVR4nO3df7BU5Z3n8feHC4HRkIjCKnIRiDETiSmQXH8kGmF1J4JlgSYZg2slmHFzR9CqpGYSo2PVxJ0sVf5K3HVrIIWrK44k6piIxCIhapQkldUIDhLAOKKBcCmQHxolGh3B7/7Rz8X20n1v//5xz+dV1XXPec5z+nz73O5vn37Oc56jiMDMzAa/Ic0OwMzMGsMJ38wsI5zwzcwywgnfzCwjnPDNzDJiaLMDABg9enRMnDix2WGYmbWVtWvX7omIMaXWb4mEP3HiRNasWdPsMMzM2oqkreXUd5OOmVlGOOGbmWWEE76ZWUa0RBt+IW+//TY9PT28+eabzQ6lrkaMGEFnZyfDhg1rdihmNsi1bMLv6elh5MiRTJw4EUnNDqcuIoK9e/fS09PDpEmTmh2OmQ1yLduk8+abb3LUUUcN2mQPIImjjjpq0P+KMbPW0LIJHxjUyb5XFl6jmbWGlk74ZmZWO074Zbjuuuu4+eabiy5fvnw5mzZtamBEZmalc8KvISd8s/qaMWMGM2bMaHYYbWvQJPxly5YxceJEhgwZwsSJE1m2bFlNnnfhwoV85CMf4cwzz+S5554D4LbbbuOUU05hypQpfO5zn+ONN97g17/+NStWrOAb3/gGU6dO5YUXXihYz8ysWQZFwl+2bBnd3d1s3bqViGDr1q10d3dXnfTXrl3LPffcw7p161i5ciVPPfUUAJ/97Gd56qmneOaZZzjxxBO5/fbb+dSnPsXs2bO56aabWLduHccff3zBemZmzTIoEv611157yNHzG2+8wbXXXlvV8/7yl7/kwgsv5LDDDuMDH/gAs2fPBmDDhg18+tOf5uMf/zjLli1j48aNBdcvtZ6ZWSO07IVX5fjDH/5QVnm1Lr30UpYvX86UKVO48847efzxx6uqZ2bWCCUf4UvqkPRvkh5K85MkPSlps6R7Jb0vlQ9P85vT8ol1iv2g4447rqzyUp111lksX76cP//5z+zbt48f//jHAOzbt4+xY8fy9ttvv6fZaOTIkezbt+/gfLF6Zlnlk67NVU6TzleBZ/PmbwBuiYgPA68Al6Xyy4BXUvktqV5dLVy4kMMOO+w9ZYcddhgLFy6s6nmnTZvGF77wBaZMmcKsWbM45ZRTAPj2t7/NaaedxhlnnMFHP/rRg/Xnzp3LTTfdxMknn8wLL7xQtJ6ZWVNExIAPoBN4FDgbeAgQsAcYmpZ/EliVplcBn0zTQ1M99ff8n/jEJ6KvTZs2HVLWn7vvvjsmTJgQkmLChAlx9913l7V+M5X7Ws3a1fTp02P69OkVrTt//vwAAoiOjo6YP39+bYNrQ8CaKCGH9z5KbcP/n8BVwMg0fxTwx4jYn+Z7gHFpehywLX2Z7Jf0aqq/p8zvorJccsklXHLJJfXchJk1yNjO49i5fVvR5QcOHGDx4sUsXrz4YNkx48azo6c+5+0GiwETvqTzgV0RsVbSjFptWFI30A3Vt7Wb2eCyc/s2JnzzoYPzW2+cDfHOoRU1hAlXrcjVueH8RoXXtko5wj8DmC3pPGAE8AHgfwFHSBqajvI7ge2p/nZgPNAjaSjwQWBv3yeNiCXAEoCurq6o9oWYWWvp7yi90KCB/R6hF0r2/ZVbQQMm/Ii4BrgGIB3hfz0iLpH0r8DngXuAecCDaZUVaf7/peU/T21NZpYhfY/SAXZ+/2oAjvmv1x9Sv98jdA0peoRvpatmb30T+DtJm8m10fdeRno7cFQq/zvg6upCNLNaalbXyD2rFvHWtg28tW0DW2+czZ5Vi0pe9/ApM8sqt8LKuvAqIh4HHk/TLwKnFqjzJvDXNYjNzAaJPasW8fq6le8WxDsH50efu2DA9XvrHHwODeHwKTNLWtfe1Ta/h8Z2Hoekmj3Gdg58ovjWW2/lxBNPZNSoUVx/fe4nqEfENCvf68/8tKzyQkafu4Dh409i+PiTmHDVioYm+8FywVjbDK1QqD2wGqWc0V+0aBGPPPIInZ2dB8uWL1/O+eefz+TJk2sWi1k76U18ZQ0V4pOuLaFtjvAb7fLLL+fFF19k1qxZ3HLLLVx55ZUFh0A2sxIUO7nawJOug+UovRpO+EV873vf49hjj+Wxxx5j1KhRAAWHQDazgfmka2tomyYdM2tfjT7pWtNrAAYRJ3yzDFmwYAGrV68GYOjQoXR3d7NoUendI6sx+twF7N+bS6qF+uGXotT1anoNwCDiJp0y9R0C2axdLFiw4D1jz/SOR7Nggbs2ZkXbHOEfM258Tb+Fjxk3vqL15s6dy1e+8hVuvfVW7r//frfjW9tYsmRJ0fJSj/Kb+QuhWQbTa26bhN+M9rUtW7YAuTtXXXrppQCcccYZ7odvbenAgQNllUP/beHtNGJl71W+kBuIrdTzB8V+FQFtmfTdpGOWER0dHWWVw7tt4RO++VC/XSt76/Q3pHGzFLvKt5ShHfr7VdSOnPDNMqK7u7us8kO06cVT1VzlW8mvolbW0k06EVGwC9Vg4oFErVF6myB6myQ6OjrKa4+uwYiVlfbOqUoVX1QdHR0Fk3t/v4paWcse4Y8YMYK9e/cO6oQYEezdu5cRI0Y0OxTLiEWLFjF9+nSmT5/O/v37y2qHbtuLp6q4yrfqX0UtpmWP8Ds7O+np6WH37t3NDqWuRowY8Z6xesxaVbuOWHn4lJnvbcPPKx9I1b+KWkzLJvxhw4YxadKkZodhZnlqcfFUo1X7RbVo0aKDPfPKGjAuT0UDztVByyZ8M7NaaccvqnoYsBFL0ghJv5H0jKSNkv57Kr9T0u8lrUuPqalckm6VtFnSeknT6vwazMysBKUc4b8FnB0Rf5I0DPiVpJ+kZd+IiPv71J8FnJAepwGL018zM2uiUm5iHsCf0uyw9Oiv68wc4K603hOSjpA0NiJ2VB2tmVWt2e3I1jwldcuU1CFpHbALeDginkyLFqZmm1skDU9l44D8y+16Ulnf5+yWtEbSmsHeE8fM2keh26muXr2a1atXV3y71FZR0knbiDgATJV0BPCApJOAa4CdwPuAJcA3gX8qdcMRsSStR1dX1+DtbG9WRKv03CjXYD/pWe7tVAca1LGVBl8rq5dORPxR0mPAzIi4ORW/Jen/Al9P89uB/KEoO1OZmVnTNOOLqtUGXyull86YdGSPpL8A/gr4naSxqUzABcCGtMoK4Eupt87pwKtuv2+MrN6zM6uvuxSFmicGerRTE0Wra7XB10o5wh8LLJXUQe4L4r6IeEjSzyWNAQSsAy5P9VcC5wGbgTeAL9c8ajMrSbnNE5Cduz81QqsNvlZKL531wMkFys8uUj+AK6oPzcysvbXa4GstO3iamVm7a7XB15zwzZqgt+fG6tWrGTp0qO8rO0gtWrSI+fPnH5zv6Ohg/vz57dFLx8yq12o9N6y+ajH4Wq34CL+Aduz1kdUjxnZ83a3Wc8Oywwl/ECh2xNgOya8a7fq6W63nhmWHE/4gkNUjxnZ93ZXcTNysFpzwB4GsHjE2+3VX2vTXaj03LDt80nYQaLW+vo3Srq+7FrfNa9dxeAa7sZ3HsXP7tqLLcwMTvNcx48azo+cP9QzroEGb8LP0geju7n5PW3Z+eaM0Y3+3wuuuVCv13LDaafUrmwdtwq9Us0e2qyRxDrYbLZcqK6+7v6PGZh8xWntxws/Tzv2js3rEmIXXXeiocef3rwYKjwDpsXCsGJ+0zdOuvT7MzErhhJ+nFr0+2vGiLWhu3O247Xa84MtsUDbpVNoO3669PqyxatX0V2kT1J5Vi3hrW+72E1tvnM3hU2Yy+lx/4djABt0RfjVXX7p/dGWydrTbzKa/PasW8fq6le8WxDu8vm4le1a19jkmaw2DLuFX82Fs9sh21SbOxx9/vOEnLltheINGv+5Kmv7KvfNUsbtOvf7MT8sqN8s3YJOOpBHAL4Dhqf79EfEtSZOAe4CjgLXAFyPiPyQNB+4CPgHsBb4QEVvqFP8hqm2Hb1avj3btIdTfF2wrx12NSpr+anZj7HinvHKzPKUc4b8FnB0RU4CpwMx0r9obgFsi4sPAK8Blqf5lwCup/JZUr2GaOU5JNUfozWwmqCbuZg9v0AxNbfpTkY9ssXKzPAO+SyLnT2l2WHoEcDZwfypfSu5G5gBz0jxp+TkqdHVInTTrw1ht00azEme1cVf7BdvM9v9ytp3fJFPo6l7IXQBW7xuBHz5lZlnlZvlK6qWTbmC+Fvgw8M/AC8AfI2J/qtIDjEvT44BtABGxX9Kr5Jp99vR5zm6gG+C44yr/cAw0dkWvxYsXH/yg5l+JWM1VjKVsu9h2y4m9N45yt12ruCtZ/8CBA2XH3ftlU812+6rVtlvl4qfe3jgHT9xqiHvpWMlKSvgRcQCYKukI4AHgo9VuOCKWAEsAurq6otLnKdY2WuqHsZoPcv66/X3Ai9XJX/+Q3hfJ4VPPO/hhHiju/lQT90Dbfk/sBRJQ0W3fOLtw27OGMOGqFQNutxS12nYrGX3uAvbvzX0RFXqPmhVTVj/8iPijpMeATwJHSBqajvI7ge2p2nZgPNAjaSjwQXInb9tGRR8iDSmaQErRtCO3KuOGKhJQM09A+uSnZdCAn2pJY9KRPZL+Avgr4FngMeDzqdo84ME0vSLNk5b/PCIqPoJvF7VoWx197gKGjz+J4eNPYsJVKxryM72pbcLNPAHpk5+WQaUc4Y8FlqZ2/CHAfRHxkKRNwD2S/gfwb8Dtqf7twL9I2gy8DMytQ9wtp13bVpsZ9+FTZhZuxmrAl00ttu3mFGs3Ayb8iFgPnFyg/EXg1ALlbwJ/XZPo2ky7tq02K+5mftm06xd0r3Z6f1nrGJRj6YA/EI1W6f5u5pdku35Bm1Vq0Cb8duXEY2b14jNUZmYZ4YRvZpYRbtKpsXZtkmlm3Fndtlmj+QjfzCwjnPDNzDLCCd/MLCOc8M3MMsIJ38wsI5zwzcwywgnfzCwjnPDNzDLCCd/MLCOc8M3MMsIJ38wsI0q5xeF4SY9J2iRpo6SvpvLrJG2XtC49zstb5xpJmyU9J+ncer4AMzMrTSmDp+0H/j4inpY0Elgr6eG07JaIuDm/sqTJ5G5r+DHgWOARSR+JiAO1DNzMzMoz4BF+ROyIiKfT9D5yNzAf188qc4B7IuKtiPg9sJkCt0I0M7PGKqsNX9JEcve3fTIVXSlpvaQ7JI1KZeOAbXmr9VDgC0JSt6Q1ktbs3r27/MjNzKwsJSd8Se8Hfgh8LSJeAxYDxwNTgR3Ad8rZcEQsiYiuiOgaM2ZMOauamVkFSkr4koaRS/bLIuJHABHxUkQciIh3gNt4t9lmOzA+b/XOVGZmZk1USi8dAbcDz0bEd/PKx+ZVuxDYkKZXAHMlDZc0CTgB+E3tQjYzs0qU0kvnDOCLwG8lrUtl/wBcLGkqEMAW4G8BImKjpPuATeR6+FzhHjpmZs03YMKPiF8BKrBoZT/rLAQWVhGXmZnVmK+0NTPLCCd8M7OMcMI3M8sIJ3wzs4xwwjczywgnfDOzjHDCNzPLCCd8M7OMcMI3M8sIJ3wzs4xwwjczywgnfDOzjHDCNzPLCCd8M7OMcMI3M8sIJ3wzs4wo5RaH4yU9JmmTpI2SvprKj5T0sKTn099RqVySbpW0WdJ6SdPq/SLMzGxgpRzh7wf+PiImA6cDV0iaDFwNPBoRJwCPpnmAWeTuY3sC0A0srnnUZmZWtgETfkTsiIin0/Q+4FlgHDAHWJqqLQUuSNNzgLsi5wngiD43PDczsyYoqw1f0kTgZOBJ4OiI2JEW7QSOTtPjgG15q/Wksr7P1S1pjaQ1u3fvLjduMzMrU8kJX9L7gR8CX4uI1/KXRUQAUc6GI2JJRHRFRNeYMWPKWdXMzCpQUsKXNIxcsl8WET9KxS/1NtWkv7tS+XZgfN7qnanMzMyaqJReOgJuB56NiO/mLVoBzEvT84AH88q/lHrrnA68mtf0Y2ZmTTK0hDpnAF8EfitpXSr7B+B64D5JlwFbgYvSspXAecBm4A3gy7UM2MzMKjNgwo+IXwEqsvicAvUDuKLKuMzMrMZ8pa2ZWUY44ZuZZYQTvplZRjjhm5llhBO+mVlGOOGbmWWEE76ZWUY44ZuZZYQTvplZRjjhm5llhBO+mVlGOOGbmWWEE76ZWUY44ZuZZYQTvplZRpRyx6s7JO2StCGv7DpJ2yWtS4/z8pZdI2mzpOcknVuvwM3MrDylHOHfCcwsUH5LRExNj5UAkiYDc4GPpXUWSeqoVbBmZla5ARN+RPwCeLnE55sD3BMRb0XE78nd5vDUKuIzM7MaqaYN/0pJ61OTz6hUNg7YllenJ5UdQlK3pDWS1uzevbuKMMzMrBSVJvzFwPHAVGAH8J1ynyAilkREV0R0jRkzpsIwzMysVBUl/Ih4KSIORMQ7wG2822yzHRifV7UzlZmZWZNVlPAljc2bvRDo7cGzApgrabikScAJwG+qC9HMzGph6EAVJP0AmAGMltQDfAuYIWkqEMAW4G8BImKjpPuATcB+4IqIOFCXyM3MrCwDJvyIuLhA8e391F8ILKwmKDMzqz1faWtmlhFO+GZmGeGEb2aWEU74ZmYZ4YRvZpYRTvhmZhnhhG9mlhFO+GZmGeGEb2aWEU74ZmYZ4YRvZpYRTvhmZhnhhG9mlhFO+GZmGeGEb2aWEU74ZmYZMWDCl3SHpF2SNuSVHSnpYUnPp7+jUrkk3Spps6T1kqbVM3gzMytdKUf4dwIz+5RdDTwaEScAj6Z5gFnk7mN7AtANLK5NmGZmVq0BE35E/AJ4uU/xHGBpml4KXJBXflfkPAEc0eeG52Zm1iSVtuEfHRE70vRO4Og0PQ7YllevJ5UdQlK3pDWS1uzevbvCMMzMrFRVn7SNiACigvWWRERXRHSNGTOm2jDMzGwAlSb8l3qbatLfXal8OzA+r15nKjMzsyarNOGvAOal6XnAg3nlX0q9dU4HXs1r+jEzsyYaOlAFST8AZgCjJfUA3wKuB+6TdBmwFbgoVV8JnAdsBt4AvlyHmM3MrAIDJvyIuLjIonMK1A3gimqDMjOz2vOVtmZmGeGEb2aWEU74ZmYZ4YRvZpYRTvhmZhnhhG9mlhFO+GZmGeGEb2aWEU74ZmYZ4YRvZpYRTvhmZhnhhG9mlhFO+GZmGeGEb2aWEU74ZmYZ4YRvZpYRA94ApT+StgD7gAPA/ojoknQkcC8wEdgCXBQRr1QXppmZVasWR/j/OSKmRkRXmr8aeDQiTgAeTfNmZtZk9WjSmQMsTdNLgQvqsA0zMytTtQk/gJ9JWiupO5UdHRE70vRO4OhCK0rqlrRG0prdu3dXGYaZmQ2kqjZ84MyI2C7pPwEPS/pd/sKICElRaMWIWAIsAejq6ipYx8zMaqeqI/yI2J7+7gIeAE4FXpI0FiD93VVtkGZmVr2KE76kwyWN7J0GPgNsAFYA81K1ecCD1QZpZmbVq6ZJ52jgAUm9z/P9iPippKeA+yRdBmwFLqo+TDMzq1bFCT8iXgSmFCjfC5xTTVBmZlZ7vtLWzCwjnPDNzDLCCd/MLCOc8M3MMsIJ38wsI5zwzcwywgnfzCwjnPDNzDLCCd/MLCOc8M3MMsIJ38wsI5zwzcwywgnfzCwjnPDNzDLCCd/MLCOc8M3MMqJuCV/STEnPSdos6ep6bcfMzEpTl4QvqQP4Z2AWMBm4WNLkemzLzMxKU68j/FOBzRHxYkT8B3APMKdO2zIzsxIoImr/pNLngZkR8d/S/BeB0yLiyrw63UB3mv1L4LmaBwKjgT11eN5qOa7ytWpsrRoXtG5srRoXtG5sxeKaEBFjSn2Sim9iXq2IWAIsqec2JK2JiK56bqMSjqt8rRpbq8YFrRtbq8YFrRtbreKqV5POdmB83nxnKjMzsyapV8J/CjhB0iRJ7wPmAivqtC0zMytBXZp0ImK/pCuBVUAHcEdEbKzHtgZQ1yajKjiu8rVqbK0aF7RubK0aF7RubDWJqy4nbc3MrPX4Slszs4xwwjczy4i2T/gDDeEgabike9PyJyVNbFBc4yU9JmmTpI2SvlqgzgxJr0palx7/2KDYtkj6bdrmmgLLJenWtM/WS5rWoLj+Mm9frJP0mqSv9anTkH0m6Q5JuyRtyCs7UtLDkp5Pf0cVWXdeqvO8pHkNiu0mSb9L/68HJB1RZN1+//d1iOs6Sdvz/l/nFVm3rkOxFInt3ry4tkhaV2Tdeu6zgnmibu+1iGjbB7kTwi8AHwLeBzwDTO5TZwHwvTQ9F7i3QbGNBaal6ZHAvxeIbQbwUBP22xZgdD/LzwN+Agg4HXiySf/bneQuLGn4PgPOAqYBG/LKbgSuTtNXAzcUWO9I4MX0d1SaHtWA2D4DDE3TNxSKrZT/fR3iug74egn/634/x/WIrc/y7wD/2IR9VjBP1Ou91u5H+KUM4TAHWJqm7wfOkaR6BxYROyLi6TS9D3gWGFfv7dbIHOCuyHkCOELS2AbHcA7wQkRsbfB2AYiIXwAv9ynOfy8tBS4osOq5wMMR8XJEvAI8DMysd2wR8bOI2J9mnyB37UtDFdlnpaj7UCz9xZbywUXAD2q5zVL0kyfq8l5r94Q/DtiWN9/DoUn1YJ30gXgVOKoh0SWpGelk4MkCiz8p6RlJP5H0sQaFFMDPJK1VboiLvkrZr/U2l+IfwGbsM4CjI2JHmt4JHF2gTivsu78h9wutkIH+9/VwZWpquqNI00Sz99mngZci4vkiyxuyz/rkibq819o94bc8Se8Hfgh8LSJe67P4aXJNFlOA/w0sb1BYZ0bENHKjmV4h6awGbbckyl2sNxv41wKLm7XP3iNyv6lbrk+zpGuB/cCyIlUa/b9fDBwPTAV2kGs6aTUX0//Rfd33WX95opbvtXZP+KUM4XCwjqShwAeBvY0ITtIwcv/EZRHxo77LI+K1iPhTml4JDJM0ut5xRcT29HcX8AC5n9T5mj00xizg6Yh4qe+CZu2z5KXepq30d1eBOk3bd5IuBc4HLklJ4hAl/O9rKiJeiogDEfEOcFuR7TVznw0FPgvcW6xOvfdZkTxRl/dauyf8UoZwWAH0nr3+PPDzYh+GWkrtgrcDz0bEd4vUOab3fIKkU8n9P+r6ZSTpcEkje6fJnezb0KfaCuBLyjkdeDXv52UjFD3iasY+y5P/XpoHPFigzirgM5JGpeaLz6SyupI0E7gKmB0RbxSpU8r/vtZx5Z/7ubDI9po5FMt/AX4XET2FFtZ7n/WTJ+rzXqvHmedGPsj1KPl3cmf5r01l/0TujQ8wglzTwGbgN8CHGhTXmeR+hq0H1qXHecDlwOWpzpXARnK9Ep4APtWAuD6UtvdM2nbvPsuPS+RuYPMC8Fugq4H/z8PJJfAP5pU1fJ+R+8LZAbxNrm30MnLnfh4FngceAY5MdbuA/5O37t+k99tm4MsNim0zufbc3vdab8+0Y4GV/f3v6xzXv6T30HpySWxs37jS/CGf43rHlsrv7H1v5dVt5D4rlifq8l7z0ApmZhnR7k06ZmZWIid8M7OMcMI3M8sIJ3wzs4xwwjczywgnfDOzjHDCNzPLiP8PhAA05t/oWvAAAAAASUVORK5CYII=\n" + }, + "metadata": { + "needs_background": "light" + } + } + ], + "source": [ + "plt.bar(\n", + " np.arange(20),\n", + " m.expected_actualdata(bestfit),\n", + " alpha=1.0,\n", + " facecolor=None,\n", + " edgecolor='k',\n", + " label='fit',\n", + ")\n", + "plt.scatter(\n", + " np.arange(20),\n", + " d[: m.config.nmaindata],\n", + " alpha=1.0,\n", + " marker='o',\n", + " c='k',\n", + " label='data',\n", + " zorder=99,\n", + ")\n", + "plt.errorbar(\n", + " np.arange(20),\n", + " d[: m.config.nmaindata],\n", + " yerr=np.sqrt(d[: m.config.nmaindata]),\n", + " marker='o',\n", + " c='k',\n", + " linestyle='',\n", + ")\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s" + ] + } + ], + "nbformat": 4, + "nbformat_minor": 2, + "metadata": { + "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.7.2-final" + }, + "orig_nbformat": 2, + "kernelspec": { + "name": "python3", + "display_name": "Python 3" + } + } +} diff --git a/src/pyhf/pdf.py b/src/pyhf/pdf.py index a25b9aff71..e5fbecc7da 100644 --- a/src/pyhf/pdf.py +++ b/src/pyhf/pdf.py @@ -66,8 +66,8 @@ def _paramset_requirements_from_modelspec(spec, channel_nbins, custom_modifiers_ spec, channel_nbins ) for req in custom_modifiers_params: - for k,v in req.items(): - _paramsets_requirements.setdefault(k,[]).append(v) + for k, v in req.items(): + _paramsets_requirements.setdefault(k, []).append(v) # build up a dictionary of the parameter configurations provided by the user _paramsets_user_configs = {} @@ -211,7 +211,7 @@ class _ModelConfig(_ChannelSummaryMixin): def __init__(self, spec, **config_kwargs): super().__init__(channels=spec['channels']) _required_paramsets = _paramset_requirements_from_modelspec( - spec, self.channel_nbins,config_kwargs.pop('custom_modifiers_params') + spec, self.channel_nbins, config_kwargs.pop('custom_modifiers_params') ) poi_name = config_kwargs.pop('poi_name', 'mu') @@ -398,7 +398,9 @@ def logpdf(self, auxdata, pars): class _MainModel: """Factory class to create pdfs for the main measurement.""" - def __init__(self, config, mega_mods, nominal_rates, custom_modifiers = None,batch_size = None): + def __init__( + self, config, mega_mods, nominal_rates, custom_modifiers=None, batch_size=None + ): self.config = config self._factor_mods = [ modtype @@ -426,13 +428,12 @@ def __init__(self, config, mega_mods, nominal_rates, custom_modifiers = None,bat ) for k, c in modifiers.combined.items() } - for i,custom in enumerate(custom_modifiers): + for i, custom in enumerate(custom_modifiers): name = f'custom_mod_{str(i).zfill(3)}' self.modifiers_appliers[name] = custom if custom.op_code == 'addition': self._delta_mods.append(name) - self._precompute() events.subscribe('tensorlib_changed')(self._precompute) @@ -515,9 +516,7 @@ def expected_data(self, pars, return_by_sample=False): pars = tensorlib.astensor(pars) deltas, factors = self._modifications(pars) - allsum = tensorlib.concatenate( - deltas + [self.nominal_rates] - ) + allsum = tensorlib.concatenate(deltas + [self.nominal_rates]) nom_plus_delta = tensorlib.sum(allsum, axis=0) nom_plus_delta = tensorlib.reshape( @@ -542,7 +541,7 @@ def expected_data(self, pars, return_by_sample=False): class Model: """The main pyhf model class.""" - def __init__(self, spec, custom_modifiers = None, batch_size=None, **config_kwargs): + def __init__(self, spec, custom_modifiers=None, batch_size=None, **config_kwargs): """ Construct a HistFactory Model. @@ -563,9 +562,11 @@ def __init__(self, spec, custom_modifiers = None, batch_size=None, **config_kwar log.info(f"Validating spec against schema: {self.schema:s}") utils.validate(self.spec, self.schema, version=self.version) # build up our representation of the specification - self.config = _ModelConfig(self.spec, custom_modifiers_params = [ - c.required_parsets for c in custom_modifiers - ], **config_kwargs) + self.config = _ModelConfig( + self.spec, + custom_modifiers_params=[c.required_parsets for c in custom_modifiers], + **config_kwargs, + ) mega_mods, _nominal_rates = _nominal_and_modifiers_from_spec( self.config, self.spec