Operations on Random Variables vs Distributions vs Random Samples What is the difference between i) random variables, ii) distributions of random variables, and iii) random samples?
While trying to figure out how to average random samples from various random variables I learned that the operations summation for random variables and distributions of random variables are different. I would like to make sure I have these terms straight and I would like to make sure the operations are valid for each term.

*

*Random variables: I think of these, non-formally, as functions. They are functions that can generate infinitely long vectors of numbers. In R, I think of random variables as rnorm or rexp. Random variables are functions that spit out numbers as defined by their distribution.


*Distributions of Random Variables: In text books these are functions (PDFs). In R these would be dnorm and dexp. An example realization of a continuous distribution on a grid of evenly spaced values would be a two column set of data where one column represents an input value and the other column represents an output frequency of the associated value.


*Random Samples: This is a vector of values output from a random variable. For example, if rnorm(n=10, 5, 1) is the random variable then the resulting vector of data is a random sample. In other words, the 10 random numbers generated from the function are the random samples. This is in contrast to random variables and distributions of random variables which are represented as functions. I do not believe that random samples must be generated from a named probability distribution. They can just be a random heap of numbers.
Operations
Random Variables

For Random Variables there is a wikipedia page Algebra Of Random Variables.
Considering two random variables X and Y, the following algebraic operations are possible:

*

*Addition: $Z=X+Y=Y+X$

*Subtraction: $Z=X-Y=-Y+X$

*Multiplication: $Z=XY=YX$

*Division (mathematics)|Division: $Z=X / Y=X \cdot (1/Y)=(1/Y) \cdot X$

*Exponentiation: $Z=X^Y=e^{Y\ln(X)}$
This notation does not seem useful. This is because I know that to sum two normals I would not just add realizations like the notations seem to indicate.
For example, to obtain a random variable $C = \frac{(A+B)}{2}$ where the random variable $A$ is :$f_A(x) = \int_{y=\infty}^{x} N(y\mid 10, 5) \, dy$ and the random variable B is:$
f_B(x) = \int_{y=\infty}^{x} N(y\mid 20, 5) \, dy$. The notation does not let the user know how to treat the parameters $\mu$ or $\sigma$ in $N(y \mid \mu, \sigma)$. It turns out that the $\mu_C$ is the average of $\mu_A$ and $\mu_C$ while the $\sigma_{C}$ is $\sigma_A$ divided by $\sqrt{2}$. The resulting random variable $C$ is $f_C(x) = \int_{y=\infty}^{x} N(y\mid 15,\frac{5}{\sqrt{2}}) \, dy$. I would not say the operations on the wikipedia page were useful to help determine the resulting random variable $C$.
So it seems that the proper, technical, way to interact with random variables is to interact with them on the distribution level as shown the in next section.
Distributions of Random Variables (PDFs or PMFs)

Given that there exists some density it is possible to use Convolutions to sum the random variables. There is a list of convolutions of probability distributions.
From Why is the sum of two random variables a convolution? the accepted answer states that how one sums random variables depends on how the random variables are represented.


*

*In terms of probability mass functions
(pmf) or probability density
functions
(pdf), it is the operation of
convolution.


*In terms of moment generating functions
(mgf), it is the (elementwise) product.


*In terms of (cumulative) distribution functions
(cdf), it is an operation closely related to the convolution.  (See
the references.)


*In terms of characteristic functions
(cf) it is the (elementwise) product.


*In terms of cumulant generating functions (cgf) it is the
sum.
The first two of these are special insofar as the box might not have a
pmf, pdf, or mgf, but it always has a cdf, cf, and cgf.

Because it is rare for me to have a need to operate on named random variables I do not tend to have a need for these operations.
What I do have a need to interact with often is random samples or realizations of random variables.
Random Samples

I do not know how to perform operations such as summation, averaging, etc. on random samples.
Example 1
Consider random sample $\textbf{a}$ and $\textbf{b}$. The data points in these random samples are $a = [1,10,100]$ and $b = [100, 10, 1]$ there are a number of ways to sum.
The way np.sum() treats the problem is that it will match the index in the summation. a[0] + b[0] then a[1] + b[1] then a[2] + b[2],$a_i + b_i$, but this misses on the tails for example a[2] + b[0] = 200. The cartesian product would capture the full width of the final distribution.
I thought of a couple ways to go about this:

*

*The Cartesian product

*The bootstrap approximation

The Cartesian product gets computationally harder for larger $n$ it also can't accept different mapping structures like covariance or past errors.
As far as I know there is not a lot of guidance on working with random samples.
The Cartesian product and the bootstrap approximation was something I thought of on the fly. Does something more rigorous exist? I am basically looking for a reference for operations on random samples.
Summary of my questions

*

*Am I grasping the terms i) random variables, ii) distributions of random variables, and iii) random samples properly?

*Am I grasping operations on i) random variables, ii) distributions of random variables properly?

*How can I perform operations such as adding, averaging, etc. on random samples? To the best of my knowledge the Cartesian product and bootstrap approximations may be possible ways to go about summing and averaging many independent random samples.


Random variables and their distributions cannot be considered separately. One defines another. To find the distribution of the function of random variables (e.g. algebraic operations) you require convolutions, as you described in your second section on the distributions.
How are the algebraic operations on random variables themselves useful? This is explained in the very first sentence of the Wikipedia article

The algebra of random variables in statistics, provides rules for the symbolic manipulation of random variables, while avoiding delving too deeply into the mathematically sophisticated ideas of probability theory. [...]

So you can do the algebraic operations on a symbolic level without bothering about the "implementation details" while having guarantees that it works as you would expect.
As for random samples, what you are describing is not correct. To find the distribution of an operation on random variables, e.g. sum, you would need to consider their joint distribution. If you used a cartesian product, you would assume that each combination of the values has equal probability. Bootstrap approximation of the distribution using the three samples you mentioned would just give you probabilities for the values $\Pr(a + b = 101) = 2/3$ and $\Pr(a + b = 20) = 1/3$, that is not the "full width" of the distribution because you would just be randomly sampling the $(a_i, b_i)$ pairs. If you meant randomly permuting the combinations, again you would be assuming that the probabilities of each combination are equally likely.
Consider the simple joint distribution
$$\begin{align}
\Pr(a=0, b=0) &= 0.2 \\
\Pr(a=0, b=1) &= 0.5 \\
\Pr(a=1, b=0) &= 0.3 \\
\Pr(a=1, b=1) &= 0 \\
\end{align}$$
In such a case you cannot claim that $\Pr(a+b = 2) > 0$, as using the proposed cartesian product approach or using random permutations, would.
With random samples, you cannot directly get the "full width" of the distribution. Depending on the problem, you could approximate the distribution of the function of the samples using some statistical model. Your example considers only three samples with two unique values for the sum, so there is not much room to build any reasonable statistical model using such data.
If you had enough data, their empirical distribution approximates the true distribution. In such a case, the distribution of the empirical results of evaluating some function over the samples would approximate the distribution of the function of the random variable. But you cannot easily go out-of-sample as you suggested.
