This is a comment to John's answer, but since comments cannot have codeblocks, let me write a CW answer.
John's approach is to let $q_1,q_2,q_3,\dotsc$ be a list of proper prime powers so that $q_i - 1 \not\mid q_n-1$ for all $i < n$. Then, $$1 - \sum_i \frac{1}{(q_i-1)} + \sum_{i < j} \frac{1}{\mathrm{lcm}(q_i-1,q_j-1)} - \sum_{i < j < k} \frac{1}{\mathrm{lcm}(q_i-1,q_j-1,q_k-1)} \pm$$ approximates the natural density. Here is the SageMath code to calculate this sum:
def relevant_powers(limit):
"""lists the proper prime powers q such that q-1 is not a multiple of x-1 for previous powers"""
powers = []
for q in srange(2,limit+1):
if q.is_prime_power() and not q.is_prime() and all([(q-1)%(x-1) > 0 for x in powers]):
powers.append(q)
return powers
def approximate_density(limit):
"""approximates the natural density of simple numbers using inclusion-exclusion principle"""
sum = 0
powers = relevant_powers(limit)
for k in range(0,len(powers)+1):
for selection in Subsets(powers,k):
sum += (-1)^k * 1/lcm([x-1 for x in selection])
print(float(sum)) # print progress
return float(sum)
Then, executing approximate_density(1000000) yields 0.46211885256102014 at some point and it does not change anymore. (Of course, it does, but SageMath's precision is not enough.)
However, executing approximate_density(5000000) yields 0.462118293064357 at some point and it does not change anymore. I am not sure what this means.