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.
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$Abs
ing0
,g1
? The result no longer satisfies the implicit relation of course. $\endgroup$