In numerical weather prediction models, point thunderstorm forecasts tend to have little predictive value beyond a few hours. Thunderstorms are difficult to predict due largely to their typically small size and correspondingly limited intrinsic predictability. We present an algorithm that predicts the probability of thunderstorm occurrence by blending multiple ensemble predictions. It combines several post-processing steps: spatial neighbourhood smoothing, dressing of probability density functions, adjusting sensitivity to model output, ensemble weighting, and calibration of the output probabilities. These operators are tuned using a machine learning technique that optimizes forecast value measured by event detection and false alarm rates. An evaluation during summer 2018 over western Europe demonstrates that the method can be deployed using about a month of historical data. Post-processed thunderstorm probabilities are substantially better than raw ensemble output. Forecast ranges from 9 hours to 4 days are studied using four ensembles: a three-member lagged ensemble, a 12-member non-lagged limited area ensemble, and two global ensembles including the recently implemented ECMWF thunderstorm diagnostic. The ensembles are combined in order to produce forecasts at all ranges. In most tested configurations, the combination of two ensembles outperforms single-ensemble output. The performance of the combination is degraded if one of the ensembles used is much worse than the other. These results provide measures of thunderstorm predictability in terms of effective resolution, diurnal variability and maximum forecast horizon.