Normalizing flows are generative machine learning models which can efficiently approximate probability distributions, using only given samples of a distribution. This architecture is used to interpolate the chiral condensate obtained from QCD simulations with five degenerate quark flavors in the HISQ action. From this a model for the probability distribution of the chiral condensate as function of lattice volume, quark mass and gauge coupling is obtained. Using the model, first order and crossover regions can be classified and the boundary between these regions can be marked by a critical mass. An extension of this model to studies of phase transitions in QCD with variable number of flavors is expected to be possible.