2
$\begingroup$

Consider random variables $X$ and $Y$, their distributions are given. $Z = f_a(X, Y)$ where $f(\cdot, \cdot)$ is a deterministic, not random function $f_a: \mathbb{R}^2 \to \mathbb{R}$ depending on a deterministic real parameter $a$. For example, let $f_a(X,Y)=\sin X + e^{aXY}$.

I've read a lot about Bayesian modelling and fitting with the help of pymc3 Python module. But how can I model $Z$ in a frequentist's manner? There's no prior data, no observations, just random variables. I am to optimize $a$ over setting $\mathbb{E}Z$ to $0$ etc.

What is the most convinient way to do this? I haven't found solution neither in scipy.optimize nor scipy.stats.

Improve this question
$\endgroup$
  • $\begingroup$ Welcome to SciComp.SE! If you're asking specifically about modelling, especially frequentist vs. Bayesina, stats.stackexchange.com might be a better place. If it's about optimization, you should be more explicit about what you're trying to optimize (it's better to give the functions and objectives you're actually interested in, to avoid answers which don't end up helping you). If I understand correctly, this is not a Bayesian problem, but a stochastic programming. I don't think you'll find ready-made functions in scipy (or anywhere else), though; you'll need to implement your own. $\endgroup$ – Christian Clason Sep 14 '16 at 7:41
  • $\begingroup$ @ChristianClason I'm interested in general approach. What's interesting, Excel atRisk program provides such an optimization $\endgroup$ – Denis Korzhenkov Sep 14 '16 at 7:55
  • $\begingroup$ The more general the approach, the more work you have to put in to apply it to your problem. Also, there's a huge difference between optimization and "setting $\mathbb{E}Z$ to 0". Finally, if I read things correctly, @Risk uses genetic algorithms -- which basically boils down to "guess values for $a$ until you get it right". There's almost always smarter approaches that use problem information (for example, derivatives of the objective function), hence it's important to specify as much information as you can. $\endgroup$ – Christian Clason Sep 14 '16 at 9:20
  • $\begingroup$ What do you mean "no observations, just random variables"? What, if anything do you know about the random variables? Do you have observations of the random variables X and Y, or of Z? If you know nothing and have no observations or data, there's nothing to do, so presumably you have something, but from your question I have no idea what that something is. Can you generate observations of Z given a particular value of the parameter $a$? $\endgroup$ – Mark L. Stone Sep 14 '16 at 12:37
  • $\begingroup$ @MarkL.Stone "their distributions are given" - imagine you have e.g. pdf for $X$ and $Y$ $\endgroup$ – Denis Korzhenkov Sep 14 '16 at 12:46
1
$\begingroup$

Mathematically, the probability density function ($\operatorname{PDF}$) for $Z$ is given by the integral: $$\operatorname{PDF}(Z) = \int \delta\left(Z - f_a(X,Y)\right) \operatorname{PDF}(X,Y)\operatorname{d}X \operatorname{d}Y.$$ If the transformation $$\left[\begin{array}{c} X \\ Y \end{array}\right] \rightarrow \left[\begin{array}{c} f_a(X,Y) \\ g_a(X,Y) \end{array}\right] $$ is one to one (if not, break it down into regions where it is, and sum) then: $$\operatorname{PDF}(Z) = \int \frac{\operatorname{PDF}(X(f_a, g_a),Y(f_a, g_a))}{\left|\frac{\partial f_a}{\partial X} \frac{\partial g_a}{\partial Y} - \frac{\partial f_a}{\partial Y} \frac{\partial g_a}{\partial X} \right|}\operatorname{d}g_a. $$ Some of the simpler choices for $g_a$ include $g_a= X$ or $g_a = Y$.

As for how I would do this, given the vague setup, I would set up a Metropolis-Hastings algorithm to produce a chain of $(Z, g_a)$ pairs from the bivariate density $$\frac{\operatorname{PDF}(X(f_a, g_a),Y(f_a, g_a))}{\left|\frac{\partial f_a}{\partial X} \frac{\partial g_a}{\partial Y} - \frac{\partial f_a}{\partial Y} \frac{\partial g_a}{\partial X} \right|},$$ and just drop the auxiliary variable, $g_a$, when returning results. You may need to implement a "sample every $N^{\mathrm{th}}$ step" limitation to reduce sample correlations. If the numbers don't have to be generated on command, but can be generated as a batch, then you can use the emcee MCMC sampling package.

Improve this answer
$\endgroup$
  • $\begingroup$ Am I right you supply to calculate $\mathbb{E}Z$ numerically with the help of samle's mean? And how should I compute such $a$ that makes $\mathbb{E}Z$ equal to $0$? $\endgroup$ – Denis Korzhenkov Sep 15 '16 at 7:29
  • $\begingroup$ For that you'll want to turn to the definition of $\mathbb{E} Z$. $$\mathbb{E} Z = \int f_a(X, Y) \operatorname{PDF}(X,Y) \operatorname{d} X \operatorname{d} Y$$ Optimizing for $\mathbb{E} Z = 0$ is the same as minimizing $(\mathbb{E} Z)^2$. With any optimization problem like this, you probably want to have a deterministic evaluation of the integral (eg scipy.integrate's dblquad), unless you want to custom build an algorithm that allows for minimization of a non-deterministic function. $\endgroup$ – Sean Lake Sep 15 '16 at 7:48
  • $\begingroup$ but what if I'd like to set, e.g., a 5% quantile of $Z$ to $0$? This problem involves inversed cdf. What's an appropriate algorithm for this case? $\endgroup$ – Denis Korzhenkov Sep 15 '16 at 7:52

Not the answer you're looking for? Browse other questions tagged or ask your own question.