# Solving Monster Equations with Mathematica

I would like to solve the following equations in Mathematica. Here are the formulas:

$u=0.5$

$\alpha=0.5$

$\epsilon_0=0.1$

$\epsilon_1=0.1$

$f_0(y)=\mbox{Gaussian density with mean -1 and variance 1}$

$f_1(y)=\mbox{Gaussian density with mean 1 and variance 1}$

$$(1-u)\left(\frac{g_1(y)}{g_0(y)}\right)^{u}-\frac{\lambda_0}{\alpha(1-\alpha)}\left(1-\alpha \Bigg(\frac{g_0(y)}{f_0(y)}\Bigg)^{\alpha-1}-(1-\alpha)\left(\frac{f_0(y)}{g_0(y)}\right)^\alpha\right)-\mu_0=0$$

$$u\left(\frac{g_1(y)}{g_0(y)}\right)^{u-1}-\frac{\lambda_1}{\alpha(1-\alpha)}\left(1-\alpha \Bigg(\frac{g_1(y)}{f_1(y)}\Bigg)^{\alpha-1}-(1-\alpha)\left(\frac{f_1(y)}{g_1(y)}\right)^\alpha\right)-\mu_1=0$$

$$\int_{-\infty}^\infty g_0(y)\mathrm{d}y=1$$

$$\int_{-\infty}^\infty g_1(y)\mathrm{d}y=1$$

$$\int_{-\infty}^\infty g_0(y)^\alpha f_0(y)^{(1-\alpha)}\mathrm{d}y=\frac{1/(\alpha (1-\alpha))- \epsilon_0}{1/(\alpha (1 - \alpha))}$$ $$\int_{-\infty}^\infty g_1(y)^\alpha f_1(y)^{(1-\alpha)}\mathrm{d}y=\frac{1/(\alpha (1-\alpha))- \epsilon_1}{1/(\alpha (1 - \alpha))}$$

In words, there are four parameters that need to be found $\lambda_0,\lambda_1,\mu_0,\mu_1$, and all other parameters are known. Whenever, the correct parameters $\lambda_0,\lambda_1,\mu_0,\mu_1$ are found, then $g_0$ and $g_1$ functions will fully be determined, which are actually the ones sought for.

The functions $g_0$ and $g_1$ are therfore parameterized by these four parameters. I tried to write a code as below, but it is wrong. I am not able to define $g_0$ and $g_1$ standalone therefore FindRoot says that it cannot read these functions. These functions are obtained as a result of the solution of the first two equations given above. My code:

Subscript[\[Epsilon], 0] = 0.1;
Subscript[\[Epsilon], 1] = 0.1;
\[Alpha] = 0.5;
u = 0.5;
f0[y_] := PDF[NormalDistribution[-1, 1], y]
f1[y_] := PDF[NormalDistribution[1, 1], y]
opts = {Method -> {Automatic, "SymbolicProcessing" -> None}, AccuracyGoal -> 8};
{g0h[y_?NumericQ, l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ], g1h[y_?NumericQ, l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ]} := FindRoot[{(1 - u) (g1[y, l0, l1, m0, m1]/g0[y, l0, l1, m0, m1])^u - (l0/(\[Alpha] (1 - \[Alpha]))) (1 - \[Alpha] (g0[y, l0, l1, m0, m1]/f0[y])^(\[Alpha] - 1) - (1 - \[Alpha]) (f0[y]/ g0[y, l0, l1, m0, m1])^\[Alpha]) - m0 == 0, u*(g1[y, l0, l1, m0, m1]/g0[y, l0, l1, m0, m1])^(u-1) - (l1/(\[Alpha] (1 - \[Alpha]))) (1 - \[Alpha] (g1[y, l0, l1, m0, m1]/f1[y])^(\[Alpha]-1) - (1 - \[Alpha]) (f1[y]/g1[y, l0, l1, m0, m1])^\[Alpha]) - m1 == 0}, {{g0[y, l0, l1, m0, m1], 1}, {g0[y, l0, l1, m0, m1], 1}}, Evaluate@opts]
h0[l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := NIntegrate[(g0[y, l0, l1, m0, m1] /. g0h[y, l0, l1, m0, m1]), {y, -Infinity, Infinity}, Evaluate@opts]
h1[l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := NIntegrate[(g1[y, l0, l1, m0, m1] /. g1h[y, l0, l1, m0, m1]), {y, -Infinity, Infinity}, Evaluate@opts]
h2[l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := NIntegrate[(g0[y, l0, l1, m0, m1] /. g0h[y, l0, l1, m0, m1])^\[Alpha]* f0[y]^(1 - \[Alpha]), {y, -Infinity, Infinity}, Evaluate@opts]
h3[l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := NIntegrate[(g1[y, l0, l1, m0, m1] /. g1h[y, l0, l1, m0, m1])^\[Alpha]* f1[y]^(1 - \[Alpha]), {y, -Infinity, Infinity}, Evaluate@opts]
{l00, l11, m00, m11} = {l0, l1, m0, m1} /. FindRoot[{h0[l0, l1, m0, m1] == 1, h1[l0, l1, m0, m1] == 1, h2[l0, l1, m0, m1] == ((1/(\[Alpha] (1 - \[Alpha]))) - Subscript[\[Epsilon], 0])/(1/(\[Alpha] (1 - \[Alpha]))), h3[l0, l1, m0, m1] == ((1/(\[Alpha] (1 - \[Alpha]))) - Subscript[\[Epsilon],1])/(1/(\[Alpha] (1 - \[Alpha])))}, {{l0, .4}, {l1, 0.4}, {m0, 0.3}, {m1, 0.3}}, StepMonitor :> Print["Step to l0,l1,m0,m1 = ", {l0, l1, m0, m1}, Evaluate@opts]]


Addendum: An intuitive (exhaustive) solution would be to choose $\lambda_0,\lambda_1,\mu_0,\mu_1$ randomly and solve the given first two equations for the variables $g_0(y)$, $g_1(y)$. After doing this for every $y$, we obtain $g_0$ and $g_1$ functions. Then check if the other $4$ equations hold or not. If not then choose another $\lambda_0,\lambda_1,\mu_0,\mu_1$, and iterate the process until the latter $4$ equations hold. But of course such an approach is very inefficient, thats the reason why there is Mathematica, hopefully..

• If $\alpha=1/2$, as specified above, the equations can be simplified quite a bit. This may allow you to pin down the form of $g_0$ and $g_1$ enough to make progress. – user46676 Mar 26 '17 at 3:35
• @marmot yes, I was aware of it after I typed the question. In such a case, I am able to obtain much simpler equations and I am also able to write a code which Mathematica can run and solve. But I am basically interested in other values of $\alpha$ too. Thats why I have $6$ equations instead of $5$ as in the case of $\alpha=0.5$. – Seyhmus Güngören Mar 26 '17 at 3:38