How to Find PDF of Transformed Random Variables Numerically? Is there a way to compute and plot the transformed random variable in Python or R?
Let's say if $X_1$ and $X_2$ are independent and identically distributed random variable with PDF defined over non-negative values of $x$ as:
$$
f(x) = \frac{x^2e^{-x}}{2} 
$$
Another random variable $Y$ is defined as $Y = X_1 + X_2$.
I know that for computing the PDF of $Y$, I can define a transformed tandom variable $T = \left(X_1+X_2, X_2\right)$ and use conservation of probability and Jacobian of the transformation.
However, what I am interested in is to know if there is a way to do this transformation and compute the PDF using any existing packages/modules of R/Python.

As far as I can tell, R doesn't have such symbolic computation capabilities; though R does have some limited symbolic computation capabilities (e.g. the command D computes symbolically the derivative of a user-provided function).
Nevertheless, when the analytical computation is hard/long, or you are sceptical about your computations, you can get an approximation via the Monte Carlo (MC) approach. The idea is pretty simple. Instead of computing the analytical PDF of $Z$, draw random numbers from its distribution. You can then use these draws to approximate whatever quantity of $Z$ you are interested in.
Example.
Suppose $X_1,X_2\sim \text{Exp}(1)$, independently. Then $Z = X_1+X_2 \sim \text{Gamma}(2, 1)$. The MC approach would be to draw $N$ values (i.e. $N\times 1$ vector) from the distribution of $X_1$ and other $N$ values (i.e. $N\times 1$ vector) from the distribution of $X_2$. The element-wise sum of these vectors is another $N\times 1$ vector drawn from $\text{Gamma}(2,1)$.
Here is the R code for this.
# number of samples, the larger the better
N <- 1e4

# draw from X1
x1 <- rexp(N)

# draw from X2 (using the same N!)
x2 <- rexp(N)

# apply your funtion; here is the sum
z = x1+x2

# Distribution of Z computed(approximated) via
# Monte Carlo
hist(z, breaks = 30, probability = T, ylim=c(0,0.4))

# The TRUE distribution of Z
plot(function(x) dgamma(x, shape = 2, rate = 1), 
     n=100, xlim=c(0,12), lwd=2, add=TRUE)


Side note 1. For the MC method to be applicable, you must be able to take random draws from $X_1,X_2$. Sometimes, this is trivial (as in the example above) but sometimes it may be less trivial. This happens when $X_1$ and $X_2$ do not have known distributions (e.g. normal, gamma, Cauchy, etc.) not there are no standard methods for generating random numbers.
Side note 2. The MC method delivers a stochastic approximation. This means that, if you run the approximation again, you will likely get slightly different answers. Analytical computation doesn't have this issue, thus my advice is to use MC as last resort or as a method for double-checking your analytical calculations.
