If you're only interested in the means and variances of linear combinations of random variables (where the means, variances, and covariances exist), then consider the following.
Suppose the means and covariances are $\mu$ $\Sigma$, respectively, then the linear combination of the multivariate random variable $X$ and another known vector $\omega$ of the same dimension has mean $$\omega.\mu$$ and the variance is $\omega'.\Sigma.\omega$. Using Mathematica notation:
n = 3; μ = Table[m[i], {i, n}]; Σ = Table[cov[Min[i, j], Max[i, j]], {i, n}, {j, n}] /. cov[i_, i_] -> σ[i]^2; ω = Table[w[i], {i, n}]; (* Mean of ω.X *) mean = ω . μ (* m[1] w[1]+m[2] w[2]+m[3] w[3] *) (* Variance of ω.X *) var = ω . Σ . ω // Simplify (* 2 cov[1,2] w[1] w[2]+2 cov[1,3] w[1] w[3]+2 cov[2,3] w[2] \ w[3]+w[1]^2 σ[1]^2+w[2]^2 σ[2]^2+w[3]^2 σ[3]^2 *)
If you're interested in more complicated functions, please include an example in your question.
Addition:
Your example in a comment:
$$Z=\prod _{i=1}^3 (a+b X_i+\epsilon_i)$$
To use replacement rules to determine $E(Z)$ and $Var(Z)$ one could construct the following:
z = Product[a + b x[i] + e[i], {i, 3}] // Expand (* a^3 + a^2 e[1] + a^2 e[2] + a e[1] e[2] + a^2 e[3] + a e[1] e[3] + a e[2] e[3] + e[1] e[2] e[3] + a^2 b x[1] + a b e[2] x[1] + a b e[3] x[1] + b e[2] e[3] x[1] + a^2 b x[2] + a b e[1] x[2] + a b e[3] x[2] + b e[1] e[3] x[2] + a b^2 x[1] x[2] + b^2 e[3] x[1] x[2] + a^2 b x[3] + a b e[1] x[3] + a b e[2] x[3] + b e[1] e[2] x[3] + a b^2 x[1] x[3] + b^2 e[2] x[1] x[3] + a b^2 x[2] x[3] + b^2 e[1] x[2] x[3] + b^3 x[1] x[2] x[3] *)
Now assuming that the random variables are all independent of each other the mean of $Z$ is
rules = {e[i_]^2 -> σ^2, e[i_] -> 0, x[i_]^2 -> σ[i]^2 + μ[i]^2, x[i_] -> μ[i]}; ez = z /.rules // FullSimplify (* (a + b μ[1]) (a + b μ[2]) (a + b μ[3]) *)
The variance is
vz = ((z^2 // Expand) /. rules) - ez^2 // FullSimplify (* σ^6 + b^2 σ^4 (μ[1]^2 + μ[2]^2 + μ[3]^2 + σ[1]^2 + σ[2]^2 + σ[3]^2) + b^4 σ^2 (μ[3]^2 σ[1]^2 + μ[3]^2 σ[2]^2 + σ[1]^2 σ[2]^2 + (σ[1]^2 + σ[2]^2) σ[3]^2 + μ[2]^2 (μ[3]^2 + σ[1]^2 + σ[3]^2) + μ[1]^2 (μ[2]^2 + μ[3]^2 + σ[2]^2 + σ[ 3]^2)) + a^4 (3 σ^2 + b^2 (σ[1]^2 + σ[2]^2 + σ[3]^2)) + b^6 ((μ[1]^2 + σ[1]^2) σ[2]^2 (μ[3]^2 + σ[3]^2) + μ[2]^2 (μ[3]^2 σ[1]^2 + (μ[1]^2 + σ[1]^2) σ[3]^2)) + 2 a^3 (2 b σ^2 (μ[1] + μ[2] + μ[3]) + b^3 (μ[3] (σ[1]^2 + σ[2]^2) + μ[2] (σ[1]^2 + σ[3]^2) + μ[1] (σ[2]^2 + σ[3]^2))) + 2 a (b σ^4 (μ[1] + μ[2] + μ[3]) + b^5 (μ[3] (μ[2] (μ[2] + μ[3]) σ[1]^2 + (μ[1] (μ[1] + μ[3]) + σ[1]^2) σ[2]^2) + (μ[2] (μ[1] (μ[1] + μ[2]) + σ[1]^2) + μ[1] σ[2]^2) σ[3]^2) + b^3 σ^2 (μ[2]^2 μ[3] + μ[1]^2 (μ[2] + μ[3]) + μ[3] (σ[1]^2 + σ[2]^2) + μ[2] (μ[3]^2 + σ[1]^2 + σ[3]^2) + μ[1] (μ[2]^2 + μ[3]^2 + σ[2]^2 + σ[3]^2))) + a^2 (3 σ^4 + 2 b^2 σ^2 ((μ[1] + μ[2] + μ[3])^2 + σ[1]^2 + σ[2]^2 + σ[3]^2) + b^4 (4 μ[1] μ[3] σ[2]^2 + (μ[1]^2 + σ[1]^2) σ[2]^2 + μ[3]^2 (σ[1]^2 + σ[2]^2) + (μ[1]^2 + σ[1]^2 + σ[2]^2) σ[3]^2 + μ[2]^2 (σ[1]^2 + σ[3]^2) + 4 μ[2] (μ[3] σ[1]^2 + μ[1] σ[3]^2))) *)
The rules would need to account for any covariances or higher order joint moments which you say you'll add after the expansion.
Mean[TransformedDistribution[a + b, {Distributed[a, x], Distributed[b, y]}]]$\endgroup$2*NormalDistribution[m,s]is not a legitimate command. If you want to avoidTransformedDistribution, consider usingExpectation. $\endgroup$