A remarkable discovery of NASA's Kepler mission is the wide diversity in the average densities of planets of similar mass. After gas disk dissipation, fully formed planets could interact with nearby planetesimals from a remnant planetesimal disk. These interactions would often lead to planetesimal accretion due to the relatively high ratio between the planet size and the hill radius for typical planets. We present calculations using the open-source stellar evolution toolkit mesa (Modules for Experiments in Stellar Astrophysics) modified to include the deposition of planetesimals into the H/He envelopes of sub-Neptunes (∼1-20 M⊕). We show that planetesimal accretion can alter the mass-radius isochrones for these planets. The same initial planet, as a result of the same total accreted planetesimal mass, can have up to 5% difference in mean densities approximately several gigayears after the last accretion due to the inherent stochasticity of the accretion process. During the phase of rapid accretion, these differences are more dramatic. The additional energy deposition from the accreted planetesimals increase the ratio between the planet's radius to that of the core during rapid accretion, which in turn leads to enhanced loss of atmospheric mass. As a result, the same initial planet can end up with very different envelope mass fractions. These differences manifest as differences in mean densities long after accretion stops. These effects are particularly important for planets that are initially less massive than ∼10 M⊕ and with envelope mass fractions less than ∼10%, thought to be the most common type of planets discovered by Kepler.