Event-by-event fluctuations caused by quantum mechanical fluctuations in the wave function of colliding nuclei in ultrarelativistic heavy ion collisions were recently shown to be necessary for the simultaneous description of $R_{AA}$ as well as the elliptic and triangular flow harmonics at high $p_T$ in PbPb collisions at the Large Hadron Collider. In fact, the presence of a finite triangular flow as well as cumulants of the flow harmonic distribution that differ from the mean are only possible when these event-by-event fluctuations are considered. In this paper we combine event-by-event viscous hydrodynamics and jet quenching to make predictions for high $p_T$ $R_{AA}$, $v_2\{2\}$, $v_3\{2\}$, and $v_2\{4\}$ in PbPb collisions at $\sqrt{s_{NN}}=5.02$ TeV. With an order of magnitude larger statistics we find that high $p_T$ elliptic flow does not scale linearly with the soft elliptical flow, as originally thought, but has deviations from perfectly linear scaling. A new experimental observable, which involves the difference between the ratio of harmonic flow cumulants at high and low $p_T$, is proposed to investigate the fluctuations of high $p_T$ flow harmonics and measure this nonlinear response. By varying the path length dependence of the energy loss and the viscosity of the evolving medium we find that $R_{AA}(p_T)$ and $v_2\{2\}(p_T)$ strongly depend on the choice for the path length dependence of the energy loss, which can be constrained using the new LHC run 2 data.