I'm trying to calculate a probability distribution of the output of a many-to-one function after assigning a probability distribution to its input. Here's a simplified question that summarises the issue:
Let's say $x$ are iid draws from $N(0,1)$ and I know that $y = f(x)$, where $f(x)$ is an analytic function with some turning points (the ratio of two polynomials with $+1000$ terms). What is the distribution of $y$?
I'm interested in the method of approaching such a question. In general my function $y = f(x)$ is quite complicated so my current solution is numerical: I (numerically) locate the turning points of $f$, and use these to divide the range of $f$ into a set of ranges such that $f$ is monotonic in each. I then estimate the probability distribution of $y$ conditional on $x$ being in each range (by sampling $f(x)$ on a grid of $x$, and weighting each output with the probability of that $x$). I then add these conditional distributions together to get an unconditional distribution of $y$.
This is messy and quite a lot of interpolation is required, so I feel I lose accuracy especially when $f(x)$ is near a turning point. Furthermore my code is not great at handling large numbers of turning points. I know I can try Monte-Carlo simulations by sampling $x$ and calculating $f(x)$ in each case, but it strikes me as even more overkill given that my probability distribution for $x$ and $f(x)$ are both analytic functions.
Are there any analytic approaches that would work here? I don't recall anything from my education that is relevant to this problem.
It seems to me that if $\phi$ is the probability density function of the normal law and $\psi$ is the probability density function of $f(X)$, you have $$ \psi(t) = \sum_{f(x_k) = t} \frac{1}{f'(x_k)} \phi(x_k) $$ You could then use a root solver on each interval of monotonicity of $f$ to find the $x_k$. Some polynomials root solvers find all the roots simultaneously.
Actually, there is a formula in wikipedia, you can start from there.