Complex organic molecules (COMs) are ubiquitously detected in the gas phase and thought to be mostly formed on icy grains. Nevertheless, there have not been any unambiguous detections of COMs larger than CH$_3$OH in ices reported thus far. Exploring this matter in greater detail has now become possible with the unprecedented possibilities offered by the James Webb Space Telescope (JWST) within the infrared (IR) spectral range with its very high sensitivity and spectral resolution in the critical 5$-$10 mu m range, the fingerprint region of oxygen-bearing COMs. In the JWST Observations of Young protoStars (JOYS+) program, more than 30 protostars are undergoing observation with the Medium Resolution Spectrograph (MRS) of the Mid-IR Instrument (MIRI). The goal of this study is to comprehensively explore the COMs ice signatures in one low- and one high-mass protostar: NGC 1333 IRAS 2A and IRAS 23385+6053, respectively. We performed global continuum and silicate subtractions of the MIRI-MRS spectra, followed by a local continuum subtraction in optical depth scale in the range around 6.8 and 8.6 mu m, the ice COM fingerprint region. We explored different choices for the local continuum and silicate subtraction. Next, we fit the observational data with a large sample of available IR laboratory ice spectra. We used the ENIIGMA fitting tool, a genetic algorithm-based code that not only finds the best fit between the lab data and the observations, but also performs a statistical analysis of the solutions, such as deriving the confidence intervals and quantifying fit degeneracy. We report the best fits for the spectral ranges between 6.8 and 8.6 mu m in NGC 1333 IRAS 2A and IRAS 23385+6053, originating from simple molecules and COMs, as well as negative ions. Overall, we find that ten chemical species are needed to reproduce the astronomical data. The strongest feature in this range (7.7 mu m) is dominated by CH$_4$, with contributions from SO$_2$ and OCN$^-$. Our results indicate that the 7.2 and 7.4 mu m bands are mostly dominated by HCOO$^-$. We also find statistically robust detections of COMs based on multiple bands, most notably, CH$_3$CHO, CH$_3$CH$_2$OH, and CH$_3$OCHO. We also report a likely detection of CH$_3$COOH. Based on the ice column density ratios between CH$_3$CH$_2$OH and CH$_3$CHO of NGC 1333 IRAS 2A and IRAS 23385+6053, we find compelling evidence that these COMs are formed on icy grains. Finally, the derived ice abundances for NGC 1333 IRAS 2A correlate well with those in comet 67P/GC within a factor of 5. Based on the high-quality JWST (MIRI-MRS) spectra, we conclude that COMs are present in interstellar ices, thus providing additional proof for the solid-state origin of these species in star-forming regions. In addition, the good correlation between the ice abundances in comet 67P and NGC 1333 IRAS 2A is fully in line with the idea that cometary COMs may be inherited from the early protostellar phases to a significant extent.