# Equation Solving with Mathematica

I need to solve equations using mathematica but I havent succeeded so far and I need help. Here is the mathematical formulation of the problem

$f_0(y)=\frac{1}{2\pi}e^{\frac{-(x+1)^2}{2}}$

$f_1(y)=\frac{1}{2\pi}e^{\frac{-(x-1)^2}{2}}$

$$g_0(y)=\left(e^{1+\mu_0+\lambda_0}f_0(y)^{-\lambda_0}\left(\lambda_1+\mu_1+\lambda_1\log\left(\frac{e^{1+\mu_0+\lambda_0}g_0(y)^{1+\lambda_0}f_0(y)^{-\lambda_0}}{f_1(y)}\right)\right)\right)^{-1/\lambda_0}$$

$$g_1(y)=\left(e^{1+\mu_0+\lambda_0}f_0(y)^{-\lambda_0}\left(\lambda_1+\mu_1+\lambda_1\log\left(\frac{g_1(y)}{f_1(y)}\right)\right)^{1+\lambda_0}\right)^{-1/\lambda_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)\log\left(\frac{g_0(y)}{f_0(y)}\right)\mathrm{d}y=0.1$$ $$\int_{-\infty}^\infty g_1(y)\log\left(\frac{g_1(y)}{f_1(y)}\right)\mathrm{d}y=0.1$$

The problem is to determine the density functions $g_0$ and $g_1$ given the density functions $f_0$ and $f_1$ as defined above.

There are $4$ equations and $4$ unknowns $\lambda_0,\lambda_1,\mu_0,\mu_1$. Normally these equations should be solvable with mathematica. The problem is that the density functions $g_0$ and $g_1$ are defined again in terms of $g_0$ and $g_1$, respectively. Therefore, one should first find $g_0$ and $g_1$ with FindRoot or maybe NSolve. After this one can use another FindRoot for $4$ equations for $4$ parameters.

I wrote the following code and it has difficulties with the choice of the starting points ($10^{-2}$ right now) of the first two FindRoots. Changing them results in different $g_0$ and $g_1$ for the same given $4$ parameters. Here is my code:

f0[y_] := PDF[NormalDistribution[-1, 1], y]
f1[y_] := PDF[NormalDistribution[1, 1], y]
opts = {Method -> {Automatic, "SymbolicProcessing" -> None}, AccuracyGoal -> 8};

lleq0[y_?NumericQ, l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := FindRoot[gg0[y, l0, l1, m0, m1] == (Exp[1 + m0 + l0]*
f0[y]^(-l0)*(l1 + m1 + l1*Log[(Exp[1 + m0 + l0]*gg0[y, l0, l1, m0, m1]^(1 + l0)*f0[y]^(-l0))/f1[y]]))^(-1/l0), {gg0[y, l0, l1, m0, m1], 10^-2}]

lleq1[y_?NumericQ, l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := FindRoot[gg1[y, l0, l1, m0,  m1] == (Exp[1 + m0 + l0]*
f0[y]^(-l0)*(l1 + m1 + l1*Log[gg1[y, l0, l1, m0, m1]/f1[y]])^(1 + l0))^(-1/l0), {gg1[y, l0, l1, m0, m1], 10^-2}]

g0[y_?NumericQ, l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := Abs[gg0[y, l0, l1, m0, m1] /. lleq0[y, l0, l1, m0, m1]]

g1[y_?NumericQ, l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] :=  Abs[gg1[y, l0, l1, m0, m1] /. lleq1[y, l0, l1, m0, m1]]

h0[l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := NIntegrate[g0[y, l0, l1, m0, m1], {y, -8, 8}, Evaluate@opts]
h1[l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := NIntegrate[g1[y, l0, l1, m0, m1], {y, -8, 8}, Evaluate@opts]
h2[l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := NIntegrate[g0[y, l0, l1, m0, m1]*Log[g0[y, l0, l1, m0, m1]/f0[y]], {y, -8, 8}, Evaluate@opts]
h3[l0_?NumericQ, l1_?NumericQ, m0_?NumericQ, m1_?NumericQ] := NIntegrate[g1[y, l0, l1, m0, m1]*Log[g1[y, l0, l1, m0, m1]/f1[y]], {y, -8, 8}, 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] == 0.1,  h3[l0, l1, m0, m1] == 0.1}, {{l0, 2}, {l1, 2}, {m0, 1}, {m1, 1}}, StepMonitor :> Print["Step to l0,l1,m0,m1 = ", {l0, l1, m0, m1}, Evaluate@opts]]


Note: $\lambda_0$ and $\lambda_1$ are supposed to be positive.

• Just as a remark: do you know or have you checked, that the implicit equations for $g_0$ and $g_1$ admit solutions for all $y \in (-\infty,\infty)$? I mean, if this is not the case, then the integration over all real values is impossible to performe, since the integral does not exist. Simple example: f[y_] := PDF[NormalDistribution[-1, 1], y] eq1 = g - (f[y]*Exp[g])^(1/1); N@Solve[eq1 == 0 && Element[y, Reals], g, Reals] Plot[Evaluate[g /. %], {y, -2, 2}] Resolve[ForAll[y, Exists[g, eq1 == 0]], Reals] Resolve[ForAll[y, y < -2 || y > -4/10, Exists[g, eq1 == 0]], Reals] – Mauricio Fernández Apr 20 '17 at 11:07
• why are you taking Abs in g0,g1 ? The result no longer satisfies the implicit relation of course. – george2079 Apr 20 '17 at 20:40
• @george2079 otherwise the search goes into complex numbers and it takes alot ot time. So this is the only reason, and normally I would not put them. – Seyhmus Güngören Apr 20 '17 at 20:41

These are only the first steps towards an answer. If you know g0 and g1 for one specific y, you know them for all y. That is, you want Mathematica to solve

Simplify[myg0 (Exp[1 + m0 + l0]*
myf0^(-l0)*(l1 + m1 +
l1*Log[(Exp[1 + m0 + l0]*myg0^(1 + l0)*myf0^(-l0))/myf1]))^(1/
l0)]==1


Mathematica claims it cannot do that, but if you do a simple variable transformation z = myg0^l0, it works

 Solve[Simplify[z (Exp[1 + m0 + l0]*
myf0^(-l0)*(l1 + m1 +
l1*Log[(Exp[1 + m0 + l0]*z^((1 + l0)/l0)*myf0^(-l0))/
myf1]))] == 1, z]


This yields

z= (E^(-1 - l0 - m0)*l0*myf0^l0)/((1 + l0)*l1*ProductLog[(E^(-1 - m0/(1 + l0))*l0*
myf0^l0)/((1 + l0)*l1*
(E^(-1 - m1/l1)*myf0^l0*myf1)^
(l0/(1 + l0)))])


which tells you that

 g0[y]=((E^(-1 - l0 - m0)*l0*f0[y]^l0)/((1 + l0)*l1*ProductLog[(E^(-1 - m0/(1 + l0))*l0*
f0[y]^l0)/((1 + l0)*l1*
(E^(-1 - m1/l1)*f0[y]^l0*f1[y])^
(l0/(1 + l0)))]))^(1/l0)


g1 can be obtained analogously. The constants can then be fixed by your integral conditions.

• I think there is a problem with the change of variables here. if $z=mygo^{l_0}$, then $mygo=z^{1/l_0}$. If you insert this into the equation, you wont get the one that you found. – Seyhmus Güngören Apr 19 '17 at 16:14
• 'mygo' is definitely a function but in your expression $mygo(a,b,c...)^{1/l_0}$, it is taken as function composition.Actually, it is just multiplication by this function, namely $mygo*(a,b,c...)^{1/l_0}$ . Therefore, as I said, the simplification is unfortunately wrong. – Seyhmus Güngören Apr 19 '17 at 16:21
• @SeyhmusGüngören Sorry, I do not understand your concerns. – marmot Apr 19 '17 at 18:32
• the answer is simply wrong – Seyhmus Güngören Apr 19 '17 at 18:33
• @SeyhmusGüngören I disagree. Please insert some random values for the constants to verify that g0 solves the third equation in your question. – marmot Apr 19 '17 at 21:48

If we assume that g0 = (h0^(-1/λ0)) f0, then we can solve for h0. We find that

coeffs = {b -> -2*E^(1 + λ0 + μ0)*λ1,
c -> -(E^(1 + λ0 + μ0)*λ1) -
(E^(1 + λ0 + μ0)*λ1)/λ0,
a -> 2*E^(1 + λ0 + μ0)*λ1 +
E^(1 + λ0 + μ0)*λ0*λ1 +
E^(1 + λ0 + μ0)*λ1*μ0 +
E^(1 + λ0 + μ0)*μ1}

h0 = -(c*ProductLog[-(E^(-(a/c) - (b*y)/c)/c)]) /. coeffs


Now that we know what g0 looks like we are confronted with the reality that MMA may not be able to integrate the function, it may not be normalized the way we want, and may not even be normalizable without a suitable weighting function.

Here is code that gives an expression for g0. We start off by defining the left hand side and the right hand side of the equation for g0

ClearAll["Global*"];

f0 = PDF[NormalDistribution[-1, 1], y];
f1 = PDF[NormalDistribution[1, 1], y];
g0 = h0^(-1/λ0) f0 ;
lhs = g0;
rhs = (Exp[1 + μ0 + λ0] f0^-λ0 (λ1 + μ1 +
λ1 Log[ Exp[1 + μ0 + λ0] g0^(1 + λ0)
f0^-λ0/f1]))^(-1/λ0);


Next, we manipulate the RHS to expand the logs of products, the logs of powers, etc.

rhs = rhs //.
Log[Times[ξ_, ζ_]] -> Log[ξ] + Log[ζ];
rhs = rhs //. Log[Power[ξ_, ζ_]] -> ζ Log[ξ];
rhs = rhs //.
Log[Times[ξ_, ζ_]] -> Log[ξ] + Log[ζ];
rhs = rhs /. Log[Exp[ξ_]] -> ξ;
rhs = rhs //.
Power[Times[ξ_, ζ_], η_] ->
Power[ξ, η] Power[ζ, η];
rhs = rhs //. Power[Exp[ξ_], η_] -> Exp[η  ξ] //
PowerExpand;


Now we multiply and exponentiate both sides to get an expression for h0 in terms of Log[h0]

{lhs, rhs} = Thread[Times[{lhs, rhs}, Sqrt[2 π]]] // Simplify;
{lhs, rhs} = Thread[Times[{lhs, rhs}, Exp[(1 + y)^2/2]]] // Simplify;
{lhs, rhs} = Thread[Power[{lhs, rhs}, -λ0]] // PowerExpand;
lhs == rhs


We note that the RHS is of the form a + b y + c Log[h0]. We find the coefficients a, b, c by

coeffs = CoefficientRules[rhs, {y, Log[h0]}];
coeffs = coeffs /. {0, 0} -> a /. {1, 0} -> b /. {0, 1} -> c;
coeffs // TableForm
rhs == a + b y + c Log[h0] /. coeffs // Simplify


Finally, solve for h0 in terms of a, b, c, y.

soln = First@Solve[lhs == a + b y + c Log[h0], h0]

g0 /. soln /. coeffs // Simplify

• Sorry, did you realize that your g0 coincides with the one I found? BTW, using the strategy I outlined yields g1[y_,l0_,l1_,m0_,m1_]=(l0/((1 + l0)*l1* ((E^(1 + l0 + m0)*(2*Pi)^(l0/2))/ (E^(-(1 + y)^2/2))^l0)^ (1 + l0)^(-1)*ProductLog[ (E^((l0*(l1 + m1 - l1*Log[1/(E^((-1 + y)^2/ 2)*Sqrt[2*Pi])]))/ ((1 + l0)*l1))*l0)/ ((1 + l0)*l1* ((E^(1 + l0 + m0)*(2*Pi)^ (l0/2))/(E^(-(1 + y)^2/ 2))^l0)^(1 + l0)^ (-1))]))^(1 + l0^(-1))` – marmot Apr 21 '17 at 15:17
• @marmot Thanks for pointing out that the solutions agree. I had not realized that before. And thanks for giving an expression for $g_1$. After plotting $g_0$ for some values of the parameters, I am thinking there may be ways to bound the solution for the parameters. – LouisB Apr 21 '17 at 19:01
• Thanks for checking. So, according to Seyhmus Güngören, your solution must be wrong ;-) (I'm kidding, both solutions are correct, I believe.) I also have problems in doing the numerical integration, but after Seyhmus Güngören's reaction my motivation to fix this dropped. – marmot Apr 21 '17 at 19:07