how to obtain the distribution of random variables? If we have two independent random variables modeled as triangular distribution centered around a mean in the interval (mean - a, mean + a). 'a' is finite value. Then discretized them at fixed sub-intervals. what will be the distribution of their max?
 A: Your question can be answered by considering the question: "in what ways can the $k$ be the maximum of these two random variables?"
Case 1: $X < k$ and $Y = k$. This happens with probability $P(X<k)P(Y=k)$, since the random variables are independent. This is indeed the first summand in $(1)$.
Case 2: $Y <k$ and $X=k$. Similarly, this happens with probability $P(Y<k)P(X=k)$. This is the second summand.
Case 3: It can also be the case (since the transformed distributions are discrete) that $Y=k$ and $X=k$. This happens with probability $P(Y=k)P(X=k)$. This is where the third summand comes into play, and there you have formula $(1)$.
In case the variables had been continuous, then the relevant formula would have been (with appropriate modifications) $(2)$, as case three happens with probability measure $0$.
A: Comment:
I think you have a good start on this (+1), but it seems some of the counts in your PDF for $Z$ are
a little off. I can't replicate them in a simple
simulation in R.
set.seed(705) # for reproducibility
p = c(1,2,3,4,3,2,1)  # proportions 
x = sample(1:7, 10^7, rep=T, prob=p)
y = sample(2:8, 10^7, rep=T, prob=p)
z = pmax(x,y)

head(cbind(x, y, z)) # first 6 of ten million
     x y z
[1,] 6 5 6
[2,] 6 4 6
[3,] 6 5 6
[4,] 7 3 7
[5,] 1 4 4
[6,] 2 4 4

round(table(z)/10^7, 4)  # tabulated probabilities
z
     2      3      4      5      6      7      8 
0.0117 0.0587 0.1642 0.2734 0.2539 0.1755 0.0624 

Your PDF is as follows:
pdf = c(7,27,42,78,79,50,16)
pdf = pdf/sum(pdf)
sum(pdf)
[1] 1
round(pdf, 4)
[1] 0.0234 0.0903 0.1405 0.2609 0.2642 0.1672 0.0535

Graph: Simulated values are heights of bars. Your PDF values are shown by (centers of) red circles.
plot(table(z)/10^7, ylab="Probability", xlab="z")
 abline(h = 0, col="green2")
 points(2:8, pdf, col="red")


Ten million iterations should give three place accuracy,
which is better than the resolution of the graph. So clear discrepancies in the graph probably represent actual
errors.
It seems especially easy to compute the exact probability $P(Z = 2) = 3/256 =   
0.0117,$ which agrees with the simulation, but not with
your $0.0234.$
I don't know if this is a source of numerical error
because I haven't checked all your arithmetic, but
it seems to me your formula
$P(X<k) = \sum_{t = -\infty}^k P(X=t)$ is not correct.
I hope some of this is helpful.
