Compute P(X=Y) when X and Y are i.i.d. geometrically distributed Let $X$ and $Y$ be independent random variables each having the same geometric distribution. We need to find $P(X=Y)$.
I tried as follows :
Let us say $X$~$Geo(p)$ and $Y$~$Geo(\theta)$.
$P(X=Y) = \sum_{y} P(X=Y|Y=y)P(Y=y) =>\sum_{y} P(X=y)P(Y=y)$.
=> $\sum_{y} p(1-p)^{y-1} \theta(1- \theta)^{y-1}$ => $\sum_{y} (p \theta)[(1-p)(1- \theta)]^{y-1} $  => $ \dfrac{p \theta}{1-(1-p)(1- \theta)}$
Is this correct ? (The independence condition is never used ). 
 A: First, there are two styles of geometric random variable. Yours
counts the trials until the first Success, the other counts
Failures until the first Success. 
Second, 'iid' means independent and identically distributed, so
that you need to have $p = \theta.$ In case $p = \theta = 1/3,$
your answer computes to $1/(9-4) = 1/5 = 0.200.$
Finally, from the discussion in the Comments, I am not entirely
certain you believe you are on the right track. So I will show you
results of a simulation in R statistical software with a million
realizations of $X$ and $Y$ with $p = 1/3.$ 
m = 10^6     # iterations
x = rgeom(m, 1/3)+1;  y = rgeom(m, 1/3)+1  # independent counts of 'trials'
mean(x == y)
## 0.200079  # Consistent with P(X = Y) = 1/5

Notes: (1) Because R uses the 'other'
version of the geometric random variable, I add $1$ to each
realization in order to match the model you are using.
(2) The comparison vector x == y has a million entries, each
of then either TRUE (for equality) or FALSE. The mean of
this vector is its proportion of TRUEs. (3) With a million
iterations, the answer should be correct to three decimal places.
So the answer is consistent with your derivation.
