1
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • $\begingroup$ 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] $\endgroup$ Commented Apr 20, 2017 at 11:07
  • $\begingroup$ why are you taking Abs in g0,g1 ? The result no longer satisfies the implicit relation of course. $\endgroup$ Commented Apr 20, 2017 at 20:40
  • $\begingroup$ @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. $\endgroup$ Commented Apr 20, 2017 at 20:41

2 Answers 2

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Commented Apr 19, 2017 at 16:14
  • $\begingroup$ '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. $\endgroup$ Commented Apr 19, 2017 at 16:21
  • $\begingroup$ @SeyhmusGüngören Sorry, I do not understand your concerns. $\endgroup$ Commented Apr 19, 2017 at 18:32
  • $\begingroup$ the answer is simply wrong $\endgroup$ Commented Apr 19, 2017 at 18:33
  • $\begingroup$ @SeyhmusGüngören I disagree. Please insert some random values for the constants to verify that g0 solves the third equation in your question. $\endgroup$ Commented Apr 19, 2017 at 21:48
0
$\begingroup$

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
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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)) $\endgroup$ Commented Apr 21, 2017 at 15:17
  • $\begingroup$ @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. $\endgroup$ Commented Apr 21, 2017 at 19:01
  • $\begingroup$ 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. $\endgroup$ Commented Apr 21, 2017 at 19:07

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.