Why is the probability of getting at least one 1 after 3 rolls of a three-sided dice not 1 - (2/3)**3? I've been told that the probability of getting at least one 1 after three rolls of a three-sided dice is $1 - (2/3)^3$. But if I try to do the math it doesn't seem true.
When I compute the probability of the event above, the sample space consists of the triples $(x, y, z)$ where the first component is the result of the first roll, the second component the result of the second roll and the third component the result of the third roll.
Now, how many of these triple exists? Maybe $10$ via $$3!/3! + 3 * 2 + 3$$ or using Python:
>>> from itertools import combinations_with_replacement
>>> list(combinations_with_replacement((1,2,3), r=3))
[(1, 1, 1), (1, 1, 2), (1, 1, 3), (1, 2, 2), (1, 2, 3), (1, 3, 3), (2, 2, 2), (2, 2, 3), (2, 3, 3), (3, 3, 3)]
>>> len(list(combinations_with_replacement((1,2,3), r=3)))
10

Of them, I'm interested in the ones where 3 appears at least once. Again, I can compute that number via $$2!/2! + 2 + 2 + 1$$ or using Python:
>>> list(filter(lambda p: 3 in p, combinations_with_replacement((1,2,3), r=3)))
[(1, 1, 3), (1, 2, 3), (1, 3, 3), (2, 2, 3), (2, 3, 3), (3, 3, 3)]
>>> len(list(filter(lambda p: 3 in p, combinations_with_replacement((1,2,3), r=3))))
6

The number I get for the overall probability is $0.6$, while I would expect to get $$1 - (2/3)^3 \approx 0.70$$
Why is that? What's wrong with my reasoning?

Your code does not count the total number of number of possible throws using the three-sided die. Rather, it counts only the possible values in those throws (i.e. ignoring the order in which they are thrown).
For example, we would consider $(1, 1, 3)$, $(1, 3, 1)$ and $(3, 1, 1)$ as distinct possibilities when we roll the die three times. Your code says that these possibilities should count as one single throw.
As @lulu notes in the comment above, there are $3^3 = 27$ possible rolls of the three-sided die. The chance that a 3 is never thrown means $2^3$ possible throws. This means that there are $27 - 8 = 19$ possible throws where a 3 can appear. This is what the formula in the title of the question counts:
$$1 - \left(\frac{2}{3}\right)^3 = \frac{19}{27} \approx 0.7$$

To emphasise why we care about order, here are the two situations phrased slightly differently: 


*

*Roll a single three-sided die three times. What is the probability that you see at least one 3 among these throws?

*Hold three three-sided dice in your hand and throw them together. What is the probability that you see a 3 on the table?
In situation 1, we care about which die shows which value. In situation 2, we don't distinguish between the dice, we just look at the values . Situation 1 is answered by $1 - \left(\frac{2}{3}\right)^3$. Your Python code and written working answers situation 2.

You can use product rather than combinations_with_replacement to generate the space of possible throws (as opposed to just the space of values that are possible with each throw). I'd probably write it like this in Python:
>>> from itertools import product
>>> possible_throws = list(product([1, 2, 3], repeat=3))
>>> len(possible_throws)
27
>>> sum(1 for throw in possible_throws if 3 in throw)
19

