What are some methods to estimate analytical PDF of random variable from an intricate expression of random variables?

Question

Suppose, I have a random variable $y$ distributed as t- distribution ($\mu$) and a random variable $x$ distributed as gamma distribution ($\alpha,\beta$), and a variable $\theta$ distributed as uniform random variable over $[-\pi,\pi]$, $w\sim\mathcal{N}(0,\sigma^2)$ is Gaussian ditributed. I need to find distribution of $h$ given as $$h=x^2|{cos}^{-1}\theta|y+y^5{|{sin}^{-1}\theta|}^2x+\frac{x}{y}+w$$ I have tried to solve for the PDF using conventional random variable transformation but the expression turns out to be so complicated that it can not be solved in analytical form.

I do not want an exact answer of this question, as I am more interested in knowing methods to use when random variable transformation are that much intricate. Numerically, may be kernel methods are used. However, what if analytical expressions are needed, approximate they may be , e.g., bounds over PDFs or some solution using perturbations? Can some general approximations to find the analytical solutions be applied? Is there a reference for the same?

Assuming all the random variables are independent.

Answer 1

One good option: kernel density estimation (KDE)

A reasonable option that is easily scalable to highly complicated cases is to generate data from the exact distribution using the stipulated transforming function and then use kernel density estimation (KDE) to get an estimated density for the transformed random variable. The main advantage of this method is that it does not require you to approximate the transforming function and so it is effective against a wide range of complicated transforms that would otherwise be intractable.

---

Computational Implementation: Here is an example of this method implemented in R using the KDE facilities in the utilities package. This generates a set of approximating probability functions (including the density function) that can be called as separate functions once loaded into the global environment. Here is an example of implementation of the present case (taking some specified values for the parameters in your problem).$^\dagger$

#Set parameters
N <- 10^4
MU    <- 2
ALPHA <- 3
BETA  <- 2
SIGMA <- 2

#Generate data from exact distribution
Y <- rt(N, df= MU)
X <- rgamma(N, shape = ALPHA, rate = BETA)
T <- runif(N, min = -pi, max = pi)
W <- rnorm(N, mean = 0, sd = SIGMA)
H <- (X^2)*abs(acos(T/pi))*Y + (Y^5)*(asin(T/pi)^2)*X + X/Y + W

#Generate kernel density estimator
KDE_H <- utilities::KDE(H, to.environment = TRUE)
KDE_H

  Kernel Density Estimator (KDE) 
 
  Computed from 10000 data points in the input 'H'
  Estimated bandwidth = 1.770201 
  Input degrees-of-freedom = Inf 
 
  Probability functions for the KDE are the following: 
 
      Density function:                   dkde * 
      Distribution function:              pkde * 
      Quantile function:                  qkde * 
      Random generation function:         rkde * 
 
  * This function is presently loaded in the global environment

You can now compute the approximating density, CDF or quantile function for the random varable $H$ at any specified input points. Below we compute the CDF over a range of points and plot the resulting approximation.

#Generate plot of (approximating) distribution functin
HH <- -100:100
PP <- pkde(HH)
plot(HH, PP, ylim = c(0, 1), type = 'l',
     main = 'Approximate Cumulative Distribution Function',
     xlab = 'H', ylab = 'Probability')

---

$^\dagger$ I will also assume that your use of $\sin^{-1}$ and $\cos^{-1}$ are with arguments expressed in radians to ensure that the generated input values are not outside the domain of these functions.