Summary Recent studies have indicated that huff ‘n’ puff (HNP) gas injection has the potential to recover an additional 30 to 70% oil from multifractured horizontal wells in shale reservoirs. Nonetheless, this technique is very sensitive to production constraints and is impacted by uncertainty related to measurement quality (particularly frequency and resolution) and lack of constraining data. In this paper, a Bayesian workflow is provided to optimize the HNP process under uncertainty using a Duvernay shale well as an example. Compositional simulations are conducted that incorporate a tuned pressure/volume/temperature (PVT) model and a set of measured cyclic injection/compaction pressure‐sensitive permeability data. Markov‐Chain Monte Carlo (MCMC) is used to estimate the posterior distributions of the model uncertain variables by matching the primary production data. The MCMC process is accelerated by using an accurate proxy model (kriging) that is updated using a highly adaptive sampling algorithm. Gaussian processes are then used to optimize the HNP control variables by maximizing the lower confidence interval (μ‐σ) of cumulative oil production (after 10 years) across a fixed ensemble of uncertain variables sampled from posterior distributions. The uncertain variable space includes several parameters representing reservoir and fracture properties. The posterior distributions for some parameters, such as primary fracture permeability and effective half‐length, are narrower, whereas wider distributions are obtained for other parameters. The results indicate that the impact of uncertain variables on HNP performance is nonlinear. Some uncertain variables (such as molecular diffusion) that do not show strong sensitivity during the primary production strongly impact gas injection HNP performance. The results of optimization under uncertainty confirm that the lower confidence interval of cumulative oil production can be maximized by an injection time of approximately 1.5 months, a production time of approximately 2.5 months, and very short soaking times. In addition, a maximum injection rate and a flowing bottomhole pressure around the bubblepoint are required to ensure maximum incremental recovery. Analysis of the objective function surface highlights some other sets of production constraints with competitive results. Finally, the optimal set of production constraints, in combination with an ensemble of uncertain variables, results in a median HNP cumulative oil production that is 30% greater than that for primary production. The application of a Bayesian framework for optimizing the HNP performance in a real shale reservoir is introduced for the first time. This work provides practical guidelines for the efficient application of advanced techniques for optimization under uncertainty, resulting in better decision making.