Stanhf: HistFactory models in the probabilistic programming language Stan

General
1

Hi Stan community,

You may find this stanhf code interesting - it’s a way to build Stan models from a declarative specification (called HistFactory) in high-energy physics. The code itself is at github and here is the paper arXiv:2503.22188.

The most unusual features are:

  • we start from a declarative specification (cf. the imperative Stan language)
  • there is an interface to access frequentist interence on a Stan model (e.g., confidence intervals, p-values etc), though note that these make a few assumptions about the model (which are reasonable for this category of models) and use functions from pyhf.

Kind regards,

Andrew

stanhf: HistFactory models in the probabilistic programming language Stan

Andrew Fowlie

In collider physics, experiments are often based on counting the numbers of events in bins of a histogram. We present a new way to build and analyze statistical models that describe these experiments, based on the probabilistic programming language Stan and the HistFactory specification. A command-line tool transpiles HistFactory models into Stan code and data files. Because Stan is an imperative language, it enables richer and more detailed modeling, as modeling choices defined in the HistFactory declarative specification can be tweaked and adapted in the Stan language. Stan was constructed with automatic differentiation, allowing modern computational algorithms for sampling and optimization.

2 Likes
2

Thanks for sharing! I’m curious what the motivation was for designing another input format—the paper didn’t really say, at least up front. Is there an existing stockpile of these models somewhere you want to support?

Is the plan to have some automated tool to do this or are you expecting people to edit JSON files?

The last part of the Stan User’s Guide shows how to do things like use a bootstrap for frequentist confidence interval estimation (though I hear it can be very unstable).

Stan is agnostic about whether constraint terms are part of the prior or likelihood, since it is only their product that matter.

The constraint terms go into the Jacobian accumulator, which can be turned off for doing maximum likelihood. With the latest release, we made the Jacobian accumulator available to users.

P.S. You might enjoy the sourcecodepro package in LaeX. It produces much more readable fixed-width code in LaTeX.

3

Hi Bob,

Thanks for your interest and comments. I will improve the paper :)

I’m curious what the motivation was for designing another input format—the paper didn’t really say, at least up front.

I didn’t create the HistFactory input format - that was created around 2012 (originally in XML and later implemented as a JSON schema). The motivation then was to simplify sharing & preserving statistical models in high-energy physics, and allow them to be changed (or ‘rescasted’ the in language used in high-enegy physics).

In high-energy physics, we use a data analysis tool called ROOT. ROOT is a data analysis and plotting framework in object-oriented C++, with a binary data format that is similar in some ways to HDF5. ROOT can be used by compiling a program or through an interactive C++ interpreter (yes, you read that right. The original one was called cint, the new one is called cling). Explicit memory management by a user is often required, sometimes even to just plot a histogram. Have you ever come across it?

ROOT began development in the mid 90s. It’s powerful, but it’s fair to say that it can be quite awkward and hard to use, and it’s not very portable. Thus, HistFactory was created as a simple way to fully declare a model that didn’t depend on ROOT or writing any C++ (even if the HistFactory models were read by ROOT). So it was portable, easily shareable, and should have a long-lifetime.

Is there an existing stockpile of these models somewhere you want to support?

The format has been adopted by a few experimental collaborations, so yes, there is now a stockpile of these models out there. And I don’t think I’d ever persuade them to switch to Stan. And I don’t think they should - the nice thing about HistFactory is that it’s implementation-agnostic. But a converter from HistFactory to Stan is nice.

Is the plan to have some automated tool to do this or are you expecting people to edit JSON files?

These JSON files are already out there in the wild & I expect experimental collaborations to keep making them. If you are using stanhf, and your model follows the HistFactory specs, you have the choice to edit the original JSON and convert it again to Stan, or to edit directly in Stan. However, if you want to depart from the HistFactory spec, you’ll need to edit the Stan file.

The last part of the Stan User’s Guide shows how to do things like use a bootstrap for frequentist confidence interval estimation (though I hear it can be very unstable).

Thanks! I will make it clearer that there are ways of doing these things already. In stanhf, this can be done without writing more Stan code or doing simulations by using asymptotic assumptions (e.g., Wilks’ theorem and similar theorems). This functionality is there so that users can keep doing things the standard way they are done in the field.

The constraint terms go into the Jacobian accumulator, which can be turned off for doing maximum likelihood. With the latest release, we made the Jacobian accumulator available to users.

Thanks. Let me think about this & check whether I need to edit my text.

4

That’s what I was guessing. We released Stan version 1 in August 2012, though BUGS and JAGS existed at the time.

That’s largely gone out of fashion in favor of templated code, because templates avoid run-time overhead of pointer chasing. We make some limited use of objects in Stan most notably for autodiff expression graph nodes, where we have to pay the virtualization penalty anyway. Otherwise, we use the curiously recursive template pattern (CRTP) when it’s possible to statically resolve function dispatch.

My guess is that people will have to switch to something like Stan when what they want to do goes beyond what the framework supports. At least that’s what we find with our own high-level tools like brms. Tools like brms make it really easy to do things it can do, but then almost impossible to do the next more complicated thing.

How do you test the accuracy of your assumptions and the corresponding results?

5

That’s a good question about checking the results. There are quite powerful theorems about the asymptotic distributions of test-statistics in this setting (e.g., Wilks’ theorem - see [1007.1727]).

There are some recommendations for best practises about checking the validity of the theorems [1911.10237], that mostly involve carefully checking that assumptions are fulfilled by the model, cf. simulations.

6

I really don’t want to sound sarcastic here, but as a software engineer, I like to verify that my algorithms work in practice as well as working in theory. The problem with asymptotic results is that we only ever have finite data. In addition, all kinds of things can go wrong moving from proofs involving real numbers to computations involving floating point representations of those numbers. We also want to test that the software has implemented the theory correctly, which we can’t do with theory alone.

7

That’s a perfectly fair comment.

The asymptotics are used because MC simulations would usually be prohibitively computationally demanding. There have been simulation studies about the applicability of asymptotics and the rate of convergence, but they aren’t performed regularly afaik.

Large physics collaborations do have statistical committees and interaction with statistics community on these topics, to try and make sure infeeenfes are reliable. You might find the PHYSTAT meeting and papers interesting in this regard: PHYSTAT