title: Trying to put together a teaching-example with pyhf
body: <p>I'm trying to learn more about pyhf and my understanding of what the goals are might be limited. I would love to fit my HEP data outside of ROOT, but I could be imposing expectations on pyhf which are not what the authors intended for it's use.</p>
  <p>I'd like to write myself a hello-world example, but I might just not know what I'm doing. My misunderstanding could also be gaps in my statistical knowledge.</p>
  <p>With that preface, let me explain what I'm trying to explore.</p>
  <p>I have some observed set of events for which I calculate some observable and make a binned histogram of that data. I hypothesize that there are two contributing physics processes, which I call signal and background. I generate some Monte Carlo samples for these processes and the theorized total number of events is close to, but not exactly what I observe.</p>
  <p>I would like to:</p>
  <ul>
  <li>Fit the data to this two process hypothesis</li>
  <li>Get from the fit the optimal values for the number of events for each process</li>
  <li>Get the uncertainties on these fitted values</li>
  <li>If appropriate, calculate an upper limit on the number of signal events.</li>
  </ul>
  <p>My starter code is below, where all I'm doing is an ML fit but I'm not sure where to go. I know it's not set up to do what I want, but I'm getting lost in the examples I find on RTD. I'm sure it's me, this is not a criticism of the documentation.</p>
  <pre class="lang-py prettyprint-override"><code>import pyhf
  import numpy as np
  import matplotlib.pyplot as plt
  
  nbins = 15
  
  # Generate a background and signal MC sample`
  MC_signal_events = np.random.normal(5,1.0,200)
  MC_background_events = 10*np.random.random(1000)
  
  signal_data = np.histogram(MC_signal_events,bins=nbins)[0]
  bkg_data = np.histogram(MC_background_events,bins=nbins)[0]
  
  # Generate an observed dataset with a slightly different
  # number of events
  signal_events = np.random.normal(5,1.0,180)
  background_events = 10*np.random.random(1050)
  
  observed_events = np.array(signal_events.tolist() + background_events.tolist())
  observed_sample = np.histogram(observed_events,bins=nbins)[0]
  
  
  # Plot these samples, if you like
  plt.figure(figsize=(12,4))
  plt.subplot(1,3,1)
  plt.hist(observed_events,bins=nbins,label='Observations')
  plt.legend()
  plt.subplot(1,3,2)
  plt.hist(MC_signal_events,bins=nbins,label='MC signal')
  plt.legend()
  plt.subplot(1,3,3)
  plt.hist(MC_background_events,bins=nbins,label='MC background')
  plt.legend()
  
  # Use a very naive estimate of the background
  # uncertainties
  bkg_uncerts = np.sqrt(bkg_data)
  
  print(&quot;Defining the PDF.......&quot;)
  pdf = pyhf.simplemodels.hepdata_like(signal_data=signal_data.tolist(), \
                                       bkg_data=bkg_data.tolist(), \
                                       bkg_uncerts=bkg_uncerts.tolist())
  
  print(&quot;Fit.......&quot;)
  data = pyhf.tensorlib.astensor(observed_sample.tolist() + pdf.config.auxdata)
  
  bestfit_pars, twice_nll = pyhf.infer.mle.fit(data, pdf, return_fitted_val=True)
  
  print(bestfit_pars)
  print(twice_nll)
  
  plt.show()
  </code></pre>