Trouble finding average case of a find max algorithm
Let me start with your main question &quot;where my mistake took place exactly?&quot; - in short, the whole point is in the wrong choice of the probability distribution  for $Pr(t_n=t)$. You are considering the number of changes in the maximum, i.e. how many times maxNum was updated, but you are using probabilities for the  location the maximum in the array.
Before starting analysis, let me say, that, for example, CLRS  Introduction to Algorithms, Third Edition, on page 120, calls formula you are using &quot;cumbersome&quot; for calculating and choose indicator function way. Donald Ervin Knuth in TAOCP Volume 1, Third Edition, 1997, pages 96-101 are using generating functions mechanism.
Now in more detail: for simplicity, let's consider the case when you have only 4 elements in the array and let us consider only such 4s in which all elements are different i.e. $(a, b, c, d)$, where all are different. I bring the pseudocode for this case from the book  Jeffrey J. McConnell, Analysis of Algorithms: an Active Learning Approach , 2001, page 3
maxNum = a
if b &gt; maxNum then
    maxNum = b
end if
if c &gt; maxNum then
    maxNum = c
end if
if d &gt; maxNum then
    maxNum = d
end if
return maxNum

where your loop is simply decomposed into separate lines and in 1-st line is not used negative infinity. As we are considering function $T$ the number of maxNum updatings, i.e. number of assignements, then, obviously $T$ have values 1,2,3,4. Probability distribution for $T$ can be counted directly and is
$Pr(T=1)$: $\frac{1}{4}=\frac{4}{16}$
$Pr(T=2)$: $\frac{1}{4}+\frac{1}{4}\cdot \frac{1}{2}+\frac{1}{4}\cdot \frac{1}{2}\cdot \frac{1}{2}=\frac{7}{16}$
$Pr(T=3)$: $\frac{1}{4}\cdot \frac{1}{2}+\frac{1}{4}\cdot \frac{1}{2}\cdot \frac{1}{2}+\frac{1}{4}\cdot \frac{1}{2}\cdot \frac{1}{2}=\frac{4}{16}$
$Pr(T=4)$: $\frac{1}{4}\cdot \frac{1}{2}\cdot \frac{1}{2}=\frac{1}{16}$
Anticipating the question of how these calculations were obtained, let me say that it is enough to simply consider the possible options for finding the maximum in the all fours and count their fractions.
To calculate general case, where we have $n$ elements in array it's more easy to find reccurence relation for $T(n)$:
$$T(n) = n \cdot T(n-1) + (n-1)!$$
and from here obtain expected value
$$\frac{ T(n)}{n!} = \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... + \frac{1}{n-1} + \frac{1}{n} = \sum_{i=1}^{n}\frac{1}{i} - 1 =\\
= \ln n + \gamma + \varepsilon_{n} -1 = \ln n +O(1)$$
where $\gamma$ is Euler–Mascheroni constant and $\varepsilon_{n} \sim \frac{1}{2n}, n\rightarrow \infty$.
Addition.
Now more about calculating probabilities in case of 4s $(a, b, c, d)$  where all are different:
The proportion of such fours in which the maximum is in the first place, i.e. maximum is $a$, is equal to the proportion of such fours in which the maximum is in second place , i.e. maximum is $b$, and is $\frac{1}{4}$. Same for cases with $c$ and $d$.
Such fours where maximum is $a$ and only they gives $1$ assignement for maxNum, i.e. $T=1$,  and so we obtain $Pr(T=1)$: $\frac{1}{4}=\frac{4}{16}$.
Such fours where maximum is $b$ gives $2$ assignements, $T=2$, but we can obtain $2$ assignements also from cases where $c$ or $d$ is maximum, so we can say, that $Pr(T=2)$ includes $\frac{1}{4}$ plus some more numbers i.e. we obtain first summand for $Pr(T=2)$.
Now let's take fours where maximum is $c$. This case can produce values $2, 3$ for $T$: when $b \gt a$, then $T=3$ and when $b \lt a$, then $T=2$. As proportion for both is $\frac{1}{4}\cdot \frac{1}{2}$, then we obtain second summand for $Pr(T=2)$ and first summand for $Pr(T=3)$.
Last we consider 4s where maximum is $d$. Hope it's easy to understand, that here we have $4$ parts with equal proportion $\frac{1}{4}\cdot \frac{1}{2}\cdot \frac{1}{2}$: 1) $b \gt a, c \gt b$ which gives $T=4$, 2) $b \gt a, c \lt b$ and $b \lt a, c \gt a$ which gives $T=3$, 3)$b \lt a, c \lt a $ which gives $T=2$. Correspondingly we obtain all remained summands.
At end let me say, that it's easy to write code for some virtual experiment, for example on python from random import uniform, I can share if needed, and check, that practical frequences approximates theoretical probabilities.