I ran the following calculation to estimate the natural density of a counting number being divisible by 3 or 5 as follows:
import numpy as np
import matplotlib.pyplot as plt
matches = []
N = 10**4
for n in range(N):
if n % 3 == 0 or n % 5 == 0:
matches.append(1)
else:
matches.append(0)
means = np.cumsum(matches) / np.arange(1, N+1)
plt.plot(means)
plt.xlabel('n')
plt.ylabel('$f_n$')
plt.show()
The calculated result for $f_{10^4} \approx 0.4667$.
Only as I was typing up this question did the site recommend I read this post which led me to this answer which suggests I simply calculate
$$1 - \left(1 - \frac{1}{3} \right)\left(1 - \frac{1}{5} \right)$$
which gets 0.46666666666666656. But I will confess I don't know where this formula comes from or why it works.

This is essentially a very simple form of the Inclusion Exclusion formula.
Basically, it is easier to ask when your divisibility criterion is NOT satisfied (and then take the complement).
In our case, roughly $\frac{2}{3}=1 -\frac{1}{3} $ of the numbers are not divisible by $3$ while roughly $\frac{4}{5}=1 -\frac{1}{5} $ of the numbers are not divisible by $5$.
By the product rule, to "miss" both is the same as multiplying those odds, namely: $\frac{2}{3}\frac{4}{5} = \frac{8}{15} = 0.5\bar{3}$.
The complement of that is about $\frac{7}{15} = 0.46\bar{7}$
Note: One has to be careful when using a distribution on infinite sets like $\mathbb{N}$ but this logic works well here. If you work in $\mathbb{Z}_{15}$ you can make this precise.