Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have $2$ equations in $2$ parameters which involve integrations. The limits of the integrals can not be defined analytically. They also require equation solving and eventually I couldnt incorporate everything in a nice way to solve. Here is my code

density functions

f0[y_] :=1/(E^((1 + y)^2/2)*Sqrt[2*Pi])
f1[y_] :=1/(E^((-1 + y)^2/2)*Sqrt[2*Pi])  

likelihood function

l[y]=E^(-(-1 + y)^2/2 + (1 + y)^2/2)

inverse of the likelihood function

ll[y_] :=Log[y]/2

The integral limits $yu$ and $yl$ which I have problems with, where $λ0$ and $μ0$ are some parameters.

FindRoot[l[yu] == Sqrt[(λ0*(λ0 + μ0*f1[yu]))/(λ0*(λ0 + (μ0 - 1)*f0[yu]))], {yu, 3}]
FindRoot[l[yl] == Sqrt[(λ0*(λ0 + (μ0 - 1)*f1[yl]))/(λ0*(λ0 + μ0*f0[yl]))], {yl, 3}]

Whenever $λ0$ and $μ0$ are known it is easy to find a solution for $yu$ and $yl$

%% my $2$ equations are as follows

h1[y_] :=Integrate[-2 + Sqrt[λ0/(λ0 + (-1 + μ0)*f0[y])] + Sqrt[(λ0 +(-1 + μ0)*f0[y])/λ0], {y, yu, Infinity}] +Integrate[-2 + Sqrt[λ0/(λ0 + μ0*f0[y])] + Sqrt[(λ0 + μ0*f0[y])/λ0], {y, -Infinity, yl}] + Integrate[-2 + Sqrt[(f1[y]*(λ0*f0[y] + λ0*f1[y]))/(f0[y]*(λ0*f1[y] + f0[y]*(λ0 + (-1 + μ0 + μ0)*f1[y])))] + Sqrt[(f0[y]*(λ0*f1[y] + f0[y]*(λ0 + (-1 + μ0 + μ0)*f1[y])))/(f1[y]*(λ0*f0[y] + λ0*f1[y]))], {y, yl, yu}]

h3[y_]:=Integrate[Sqrt[(λ0*f0[y]^2)/(λ0 + (-1 + μ0)*f0[y])], {y, yu, Infinity}] + Integrate[Sqrt[(λ0*f0[y]^2)/(λ0 + μ0*f0[y])], {y, -Infinity, yl}] + Integrate[Sqrt[(f0[y]*f1[y]*(λ0*f0[y] + λ0*f1[y]))/(λ0*f1[y] + f0[y]*(λ0 + (-1 + μ0 + μ0)*f1[y]))], {y, yl, yu}]

I am using Findroot to solve them as follows:

FindRoot[{h1[y] == 0.08, h3[y] == 1}, {{μ0, 0.4}, {λ0, 0.8}}, 
         StepMonitor :> Print["Step to μ0 = ", μ0]]

as one can see $h1$ and $h3$ involve $yl$ and $yu$ which are the solution of another equations. These equations can only be numerically evaluated. At this point my eventual aim is to solve the last FindRoot of my code.

In simple words if I would be Mathematica I would guess some $λ0$ and $μ0$ then I would go to the first $2$ FindRoots of my code solve them and get $yl$ and $yu$ numerically.

Then using these values I would obtain $h1$ and $h3$ and eventually I would go to last FindRoot of my code to see if $h1=0.08$ and $h3=1$ are satisfied or not.

If not I would change $λ0$ and $μ0$ until the equations are satisfied. This is of course the worst way of solving the equations but Mathematica would have a better way I suppose.

How should I modify my code such that Mathematica could start solving them? The equations are horrible and will take at any case alot of time.

Thanks for reading my problems and for any help.

share|improve this question
    
You have an undefined lambda1 and mu1 in h1 and h3 –  belisarius Nov 24 '13 at 20:50
    
@belisarius my apologies $μ1:=μ0$ and $λ1:=λ0$ I just corrected. –  Seyhmus Güngören Nov 24 '13 at 21:08

1 Answer 1

up vote 2 down vote accepted

You have too many errors to describe each one of them separately.

I'll post a working code below, but you should try to learn how to solve problems with Mathematica one step at a time. You had thrown in too much non-working code and your debugging will be always a nightmare programming like that.

f0[y_] := 1/(E^((1 + y)^2/2)*Sqrt[2*Pi])
f1[y_] := 1/(E^((-1 + y)^2/2)*Sqrt[2*Pi])
l[y_] := E^(-(-1 + y)^2/2 + (1 + y)^2/2)
opts = {Method -> {Automatic, "SymbolicProcessing" -> None}, AccuracyGoal -> 8};

yueq[λ0_?NumericQ, μ0_?NumericQ] := FindRoot[l[yu] == Sqrt[(λ0*(λ0 + μ0* f1[yu]))/
                                                           (λ0*(λ0 + (μ0 - 1)* f0[yu]))], {yu, 3}]
yleq[λ0_?NumericQ, μ0_?NumericQ] := FindRoot[l[yl] == Sqrt[(λ0*(λ0 + (μ0 - 1)* f1[yl]))/
                                                          (λ0*(λ0 + μ0*f0[yl]))], {yl, 3}]

h1[λ0_?NumericQ, μ0_?NumericQ] := 
 NIntegrate[-2 + Sqrt[λ0/(λ0 + (-1 + μ0)*f0[y])] + Sqrt[(λ0 + (-1 + μ0)*f0[y])/λ0], 
                                 {y, yu /. yueq[λ0, μ0], Infinity}, Evaluate@opts] +
  NIntegrate[-2 + Sqrt[λ0/(λ0 + μ0*f0[y])] + Sqrt[(λ0 + μ0*f0[y])/λ0], 
                                 {y, -Infinity, yl /. yleq[λ0, μ0]}, Evaluate@opts] + 
  NIntegrate[-2 + Sqrt[(f1[y]*(λ0*f0[y] + λ0*f1[y]))/(f0[y]*(λ0*f1[y] + 
                       f0[y]*(λ0 + (-1 + μ0 + μ0)*f1[y])))] + 
                  Sqrt[(f0[y]*(λ0*f1[y] + f0[y]*(λ0 + (-1 + μ0 + μ0)*f1[y])))/
                       (f1[y]*(λ0*f0[y] + λ0*f1[y]))], 
                                   {y, yl /. yleq[λ0, μ0], yu /. yueq[λ0, μ0]}, Evaluate@opts]

h3[λ0_?NumericQ, μ0_?NumericQ] := 
 NIntegrate[ Sqrt[(λ0*f0[y]^2)/(λ0 + (-1 + μ0)*f0[y])], 
                                    {y, yu /. yueq[λ0, μ0], Infinity}, Evaluate@opts] + 
  NIntegrate[Sqrt[(λ0* f0[y]^2)/(λ0 + μ0*f0[y])], 
                                    {y, -Infinity, yl /. yleq[λ0, μ0]}, Evaluate@opts] + 
  NIntegrate[Sqrt[(f0[y]* f1[y]*(λ0*f0[y] + λ0*f1[y]))/(λ0* f1[y] + f0[y]*(λ0 + (-1 + μ0 + μ0)*f1[y]))],
                                     {y, yl /. yleq[λ0, μ0], yu /. yueq[λ0, μ0]}, 
   Evaluate@opts]

sol = FindRoot[{h1[λ0, μ0] == 0.08, h3[λ0, μ0] == 1}, {{μ0, 0.4}, {λ0, 0.8}}, 
               StepMonitor :> Print["Step to μ0 = ", μ0]]

(*
 {μ0 -> 0.125521 - 6.10413*10^-22 I, λ0 -> 0.286823 + 2.85401*10^-23 I}
*)
share|improve this answer
    
that is very nice of you. I will read it tomorrow carefully and let you know what I learned. You are right. My codes are not so nice. –  Seyhmus Güngören Nov 24 '13 at 22:12
    
it is really great and so fast! What I learned is that if I will not use symbolic things then I should better turn it off and the accuracygoal is also important. I think both factors make the code fast. I guess symbolic thing has higher impact of the speed of the code? I also understood that whenever I have equations to solve and there are some other parameters which can be found analytically I can state them inside findroot with "/." I also learned Evaluate@... so to use the reasonable options. I also havent understood yet the meaning of "μ0_?" the question marking. –  Seyhmus Güngören Nov 26 '13 at 15:38
    
@SeyhmusGüngören Glad it helped. See reference.wolfram.com/mathematica/ref/PatternTest.html. –  belisarius Nov 26 '13 at 15:43
    
thanks again I checked it. For all $\mu 0$ which is numeric the expression becomes true. Is it used to restrict the search to numerical values? –  Seyhmus Güngören Nov 26 '13 at 15:51
    
@SeyhmusGüngören It's used to restraint the use of symbolic processing, which may give exact results but will cost you processing time and/or errors if the symbolic results can't be obtained –  belisarius Nov 26 '13 at 15:54

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.