Re: FindRoot within a FindRoot Problem
- To: mathgroup at smc.vnet.net
- Subject: [mg84234] Re: [mg84212] FindRoot within a FindRoot Problem
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 14 Dec 2007 15:00:16 -0500 (EST)
- References: <200712140856.DAA07770@smc.vnet.net>
stephen.a.mccormack at baesystems.com wrote:
> I need to find mu1 and sigma1 such that the two errors are zero (or at
> least minimized). The following code works but is manual and iterative
> in that I need to manually guess the values for mu1 and sigma1, check
> the errors and guess again, and so on.
>
> \[Mu]1 = 80.(* a guessed value *);
> \[Sigma]1 = 107.0 (* a guessed value *);
> \!\(TraditionalForm\`\(\(solution1 = FindRoot[{\[Mu]1 ==
> \*SuperscriptBox[\(E\), \(M1 +
> \*FractionBox[
> SuperscriptBox[\(S1\), \(2\)], \(2\)]\)], \[Sigma]1 ==
> \*SqrtBox[\(
> \*SuperscriptBox[\(E\), \(
> \*SuperscriptBox[\(S1\), \(2\)] + 2*M1\)]*\((
> \*SuperscriptBox[\(E\),
> SuperscriptBox[\(S1\), \(2\)]] -
> 1)\)\)]}, {M1, .1}, {S1, .1}];\)\[IndentingNewLine]
> \(m1 = ReplaceAll[M1,
> solution1[\([\)\(1\)\(]\)]];\)\[IndentingNewLine]
> s1 = ReplaceAll[S1, solution1[\([\)\(2\)\(]\)]]\)\);
> dist1 = LogNormalDistribution[m1, s1];
> dist1pdf = PDF[dist1, x];
> ivalue2a = Integrate[dist1pdf, {x, 0, 10}]*100
> ivalue13a = 100 - Integrate[dist1pdf, {x, 0, 1000}]*100
> error1 = (ivalue2a - 2.78*10^-16)/(2.78*10^-16)*100
>
> error2 = (ivalue13a - .2)/.2*100
> {m1, s1}
>
> I tried wrapping the code within FindRoot but get error messages (see
> below). Does anyone know how to force Mathematica to solve for mu1 and
> sigma1 while setting both errors to zero ?
>
> Thanks
>
> \!\(TraditionalForm\`\(\(solution1 = FindRoot[{\[Mu]1 ==
> \*SuperscriptBox[\(E\), \(M1 +
> \*FractionBox[
> SuperscriptBox[\(S1\), \(2\)], \(2\)]\)], \[Sigma]1 ==
> \*SqrtBox[\(
> \*SuperscriptBox[\(E\), \(
> \*SuperscriptBox[\(S1\), \(2\)] + 2*M1\)]*\((
> \*SuperscriptBox[\(E\),
> SuperscriptBox[\(S1\), \(2\)]] -
> 1)\)\)]}, {M1, .1}, {S1, .1}];\)\[IndentingNewLine]
> \(m1 = ReplaceAll[M1,
> solution1[\([\)\(1\)\(]\)]];\)\[IndentingNewLine]
> s1 = ReplaceAll[S1, solution1[\([\)\(2\)\(]\)]]\)\);
> dist1 = LogNormalDistribution[m1, s1];
> dist1pdf = PDF[dist1, x];
> ivalue2a = Integrate[dist1pdf, {x, 0, 10}]*100
> ivalue13a = 100 - Integrate[dist1pdf, {x, 0, 1000}]*100
> solution2 =
> FindRoot[{0. == (ivalue2a - 2.78*10^-16)/(2.78*10^-16)*100.,
> 0. == (ivalue13a - .2)/.2*100. == 0.}, {\[Mu]1, 95.}, {\[Sigma]1,
> 107.}]
I find this quite hard to follow. In general it is vastly easier on the
eyes if you use plain text ascii instead of the human-unreadable mess
above. To do this, you might do Cell > Convert To > InputForm before you
copy/paste into email.
What you want to do is define these various entities as functions of mu
and sigma. Start with solution1[mu,sigma]. It will do a FindRoot,
returning the {m1,s1} pair appropriate for the given values of mu and
sigma. We will use this in dist1, etc.
Here is the code I used, and the resulting {mu,sigma} pair.
solution1[(mu_)?NumberQ, (sigma_)?NumberQ] :=
{m1, s1} /. FindRoot[{mu == E^(m1 + s1^2/2),
sigma == Sqrt[E^(s1^2 + 2*m1)*(E^s1^2 - 1)]}, {m1, 0.1},
{s1, 0.1}]
dist1[(mu_)?NumberQ, (sigma_)?NumberQ] :=
Module[{m, s}, {m, s} = solution1[mu, sigma];
LogNormalDistribution[m, s]];
dist1pdf[(mu_)?NumberQ, (sigma_)?NumberQ] :=
PDF[dist1[mu, sigma], x];
ivalue2a[(mu_)?NumberQ, (sigma_)?NumberQ] :=
NIntegrate[dist1pdf[mu, sigma], {x, 0, 10}]*100
ivalue13a[(mu_)?NumberQ, (sigma_)?NumberQ] :=
100 - NIntegrate[dist1pdf[mu, sigma], {x, 0, 1000}]*100
error1[(mu_)?NumberQ, (sigma_)?NumberQ] :=
((ivalue2a[mu, sigma] - 2.78/10^16)/(2.78/10^16))*100
error2[(mu_)?NumberQ, (sigma_)?NumberQ] :=
((ivalue13a[mu, sigma] - 0.2)/0.2)*100
FindMinimum[error1[mu, sigma]^2 + error2[mu, sigma]^2,
{mu, 80}, {sigma, 107}]
Out[16]= {3.24884*10^-23, {mu -> 342.766, sigma -> 142.691}}
Daniel Lichtblau
Wolfram Research
- References:
- FindRoot within a FindRoot Problem
- From: stephen.a.mccormack@baesystems.com
- FindRoot within a FindRoot Problem