Equation Solving with Mathematica

Question

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.

Answer 1

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.

Answer 2

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