We present a Bayesian hierarchical model which enables a joint fit of the ultra-high-energy cosmic ray ( UHECR ) energy spectrum and arrival directions within the context of a physical model for the UHECR phenomenology . In this way , possible associations with astrophysical source populations can be assessed in a physically and statistically principled manner . The importance of including the UHECR energy data and detection effects is demonstrated through simulation studies , showing that the effective GZK horizon is significantly extended for typical reconstruction uncertainties . We also verify the ability of the model to fit and recover physical parameters from CRPropa 3 simulations . Finally , the model is used to assess the fraction of the the publicly available dataset of 231 UHECRs detected by the Pierre Auger Observatory ( PAO ) which are associated with the Fermi -LAT 2FHL catalogue , a set of starburst galaxies and Swift -BAT hard X-ray sources . We find association fractions of 9.5 ^ { +2.4 } _ { -5.9 } , 22.7 ^ { +6.6 } _ { -12.4 } and 22.8 ^ { +6.6 } _ { -8.0 } per cent for the 2FHL , starburst galaxies and Swift -BAT catalogues respectively .