I have the following system
$$ \mathbf M \left(\alpha_1,\alpha_2,...,\alpha_n \right) \vec v = \vec b $$
$$ v_{out} = \vec c^T \vec v $$
where $\mathbf M \left(\alpha_1,\alpha_2,...,\alpha_n \right)$ is a matrix where some elements are Berboulli distributed. I wish to find or approximate the distribution of $v_{out}$
$\mathbf M$ depends linearly on the $\alpha$'s and is non-singular for $\alpha$'s in the range $\left[0;1\right]$. The dimension $n$ of $\vec \alpha$ is as large as I can handle and at least $n=100$. The $\alpha$'s are not independent, but would have some decaying covariance with their neighbors (symmetric covariance matrix, concentrated near the diagonal).
I can in principle find the derivatives of any order of $v_{out}$ against $\vec \alpha$, and I thought maybe I could find the moments (or cumulants (which I don't really know much about)) using the derivatives. However, it becomes difficult to manage when going beyond the 2nd moment for large $n$.
Based on Monte-Carlo simulations it it not sufficient to assume $v_{out}$ is normal distributed, which I assume means I need at least the 3rd moment.
The system is an electronic system with random switches, that is an approximate model of the physical system.
What is the best approach to do this?