Re: using FindMinimum to minimize a function involving NIntegrate
- To: mathgroup at smc.vnet.net
- Subject: [mg15573] Re: using FindMinimum to minimize a function involving NIntegrate
- From: Robert Knapp <rknapp>
- Date: Tue, 26 Jan 1999 13:44:50 -0500 (EST)
- Organization: Wolfram Research, Inc.
- References: <77v2ue$kle@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Albert Maydeu-Olivares wrote:
>
> Hi, everyone
>
> I can't get FindMinimum to work when the function to be minimized
> involves numerical integrals.
> Any suggestions would be most welcome.
>
> I enclose a typical example: I'm trying to minimize e .e where e =
> {data1-definite integral of a univariate normal distribution,
> data2- definite integral of a univariate normal distribution,
> data3- definite integral of a bivariatye normal distribution} using
> Newton's method.
>
> Below is the code. The gradient is correct.
>
> In[1]:=
> dens[x_]=1/(E^(x^2/2)*Sqrt[2*Pi]);
> p[h_]:=Module[{x},NIntegrate[dens[x],{x,h,Infinity}]] dens2[h_,k_,x_]=
> 1/(2*E^((h^2 + k^2 - 2*h*k*x)/(2*(1 - x^2)))*Pi*Sqrt[1 - x^2]);
> pp[h_,k_,r_]:= Module[{x},
> NIntegrate[dens2[h,k,x],{x,0,r},PrecisionGoal->10]+p[h] p[k]]
>
> In[3]:=
> nvars=2;ncorr=1;
> prop={.828,.628,.567};
> t={t1,t2};
> r={l1 *l2};
> gradient={{-dens[t1],0,0,0},
> {0,-dens[t2],0,0},
> {dens[t1]*(-1 + Erf[(-(l1*l2*t1) + t2)/Sqrt[2 - 2*l1^2*l2^2]])/2,
> dens[t2]* (-1 + Erf[(t1 - l1*l2*t2)/Sqrt[2 - 2*l1^2*l2^2]])/2,
> dens2[t1,t2,l1*l2]*l2, dens2[t1,t2,l1*l2]*l1}};
>
> In[4]:=
> e := prop-{p[t[[1]]],p[t[[2]]],pp[t[[2]],t[[1]],r[[1]]]} ; FindMinimum[e
> . e, {t1,-.95},{t2,-.4},{l1,.5},{l2,.6},
> Gradient->gradient, Method->Newton]
>
> Out[4]=
> "Could not symbolically find the gradient of (e . e)
> Try giving two starting values for each variable."
>
> I believe the error message given by FindMinimum is misleading.
>
> Trying to figure out what's going on, I specify the problem using
> Integrate and then
> numerically evaluating the integrals inside FindMinimum. It did not work
> either.
> It just sit there.
>
> In[11]:=
> p2[h_]:=Integrate[dens[x],{x,h,Infinity}]
> pp2[h_,k_,r_]:=Integrate[dens2[h,k,x],{x,0,r}]+p[h] p[k]
>
> e2 := prop-{p2[t[[1]]],p2[t[[2]]],pp2[t[[2]],t[[1]],r[[1]]]} ;
> FindMinimum[N[e2 . e2], {t1,-.95},{t2,-.4},{l1,.5},{l2,.6},
>
> --
>
Your problem here is caused by the way the Gradient option is
interpreted. You are specifying the Hessian, not the gradient. Since
in Version 3, there is no Hessian option, FindMinimum is designed to
work by looking at the TensorRank of the value and deciding whether it
is a gradient or Hessian. (This is why the example you gave in a
another post several days ago worked even though you were specifying
the Hessian there too.) In retrospect, this was a bit silly because it
leads to confusion in cases such as yours. In the next version of
Mathematica released, there will be a Hessian option.
Since that does not help you yet, I will make a couple of suggestions.
First, using the same technique as you used to find the Hessian, find
the gradient, then use Method->QuasiNewton. Typically for machine
numbers the quasi Newton method (BFGS) is faster anyway.
Second, this is a slower method, but it seems to work. Let the method
remain the default, and specify two starting values for each variable.
This does give a result:
In[7]:=
e := prop-{p[t[[1]]],p[t[[2]]],pp[t[[2]],t[[1]],r[[1]]]} ; FindMinimum[e
.. e, {t1,-.95,-1.},{t2,-.4,-.5},{l1,.5,.4},{l2,.6,.5}]
Out[7]=
-12
{2.04392 10 , {t1 -> -0.94629, t2 -> -0.32656,
l1 -> 0.762022, l2 -> 0.592751}}