Introducing constraints in Nonlinear Fit "a la Levenberg Marquardt"
- To: mathgroup at smc.vnet.net
- Subject: [mg43501] Introducing constraints in Nonlinear Fit "a la Levenberg Marquardt"
- From: gilmar.rodriguez at nwfwmd.state.fl.us (Gilmar Rodríguez Pierluissi)
- Date: Thu, 18 Sep 2003 05:39:16 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
A paper in : http : // www.primepuzzles.net/puzzles/GoldbachPartitions.pdf
conjectures that if G(n) = Length[{(p, q) : n = p + q, p <= q}] then
the lower and upper bounding functions of G(n) can be expressed
respectively as simple exponential of the form Exp[alpha*n^beta],
where alpha>0, and 0<beta<1.
Define G(n) (using Mathematica commands) as follows :
In[1]: ClearAll@G
In[2]: G[n_Integer?Positive] := {n, Length[ Transpose@{n - #, #} &@
Select[n - Prime@Range@PrimePi[n/2],
PrimeQ]]}
Notice that the function G[n] above gives a point{n, G(n)}; i.e.
not just {G(n)}.
Now compute :
In[3]: A = Table[G[n], {n, 4, 1000, 2}];
and plot A :
In[4]: plt0 = ListPlot[A, PlotJoined -> True, PlotRange -> All,
PlotStyle -> Hue[.4]]
Since A is a large set; we build a set M containing the maximal
boundary points of A:
In[5]: SMaxima[B_] := Module[{M}, {M = {}; AppendTo[M, B[[1]]];
{Xg, Yg} = B[[1]];
Do[If[Yg < B[[i + 1]]]],
{{Xg, Yg} = B[[i + 1]],
AppendTo[M, B[[i+ 1]]]}],
{i, 1, Length[B] - 1}]}; M]
In[6]: M = SMaxima[A]
and you get:
Out[6]: M =
{{4, 1}, {10, 2}, {22, 3}, {34, 4}, {48, 5}, {60, 6},
{78, 7},{84, 8}, {90, 9},{114, 10}, {120, 12}, {168, 13},
{180, 14},{210, 19}, {300, 21}, {330, 24},{390, 27}, {420, 30},
{510, 32},{630, 41}, {780, 44}, {840, 51}, {990, 52}}
Plot M :
In[7]: plt1 = ListPlot[M, PlotJoined -> True, PlotRange -> All,
PlotStyle -> Hue[.6]]
Compare A and M :
In[8]: Show[plt0, plt1]
Pick the Maximum point in M given by
{M[[Length[M]]][[1]], M[[Length[M]]][[2]]} :
In[9]: Off[General::"spell1"]
In[10]: Data = {{M[[Length[M]]][[1]], M[[Length[M]]][[2]]}}
Get:
Out[10]: Data = {{990, 52}}
A paper in: http : // www.maths.ex.ac.uk/~mwatkins/zeta/wolfgas.htm
gives a clue about what a bounding function for A might be:
In[11]: << Statistics`NonlinearFit`
In[12]: Clear[BF, alpha]
In[13]: BF[N_] = NonlinearFit[Data, alpha/(Exp[2/Log[N]] - 1),
{N}, {alpha}]
Out[13]: BF[N_] =17.4909/(-1 + [Exp[2 / Log[N])
Build the bounding set of M with points whose x-coordinates are
the same, but whose y-coordinates are approximated using the function BF[N]:
In[14]: BoundSet = Table[{M[[i]][[1]], BF[M[[i]][[1]]]}, {i, 1, Length[M]}]
Get:
Out[14]:
BoundSet =
{{4, 5.41163}, {10, 12.6421}, {22, 19.2236}, {34, 22.9164},
{48, 25.8596},{60, 27.7706}, {78, 30.0226}, {84, 30.6597},
{90, 31.2531}, {114, 33.2883},{120, 33.7304}, {168, 36.6333},
{180, 37.2292},{210,38.5612}, {300, 41.6466},{330, 42.4718},
{390, 43.9188}, {420, 44.5609}, {510, 46.244}, {630, 48.0767},
{780, 49.9301}, {840, 50.5734}, {990, 52.}}
Plot BoundSet:
In[15]: plt2 = ListPlot[BoundSet, PlotJoined -> True,
PlotRange -> All, PlotStyle -> Hue[.4]]
Plot A, M, and BoundSet together:
In[16]: Show[plt0, plt1, plt2]
On the other hand; when I attempt to find a model as the one originally
suggested in :
http : // www.primepuzzles.net/puzzles/GoldbachPartitions.pdf
I do:
In[17]: Data2 = {M[[1]], M[[Length[M]]]}
Get:
Out[17]: Data2 = {{4, 1}, {990, 52}}
In[18]: Clear[BF2, alpha, beta]
In[19]: $IterationLimit = Infinity
In[20]: BF2[N_] = NonlinearFit[Data2, alpha*Exp[N^beta], {N},
{alpha, beta}]
I get a message: "FindMinimum::fmlim:The minimum could not be
bracketed in 30 iterations."
Out[20]: BF2{N_] = 0.0187577*Exp[N^0.995686]
To complete our analysis we want to finish with:
In[21]: BoundSet2 = Table[{M[[i]][[1]], BF2[M[[i]][[1]]]}, {i, 1, Length[M]}]
In[22]: plt3 = ListPlot[BoundSet2, PlotJoined -> True,
PlotRange -> All, PlotStyle -> Hue[.4]]
and finally:
In[23]: Show[plt0, plt1, plt3]
My questions are:
(1.) How can I improve the In[20] command:
BF2[N_] = NonlinearFit[Data2, alpha*Exp[N^beta], {N}, {alpha, beta}]
so that I can get a better convergence, and get a better answer than Out[20]?
(2.) Is there any way to enter the constraints alpha>0, 0<beta<1 in
NonlinearFit to get a better answer? Or to give an initial "kick" to the
parameters alpha, and beta ( "a la Levenberg-Marquardt") to get a better
result for BF2[N]?
Thank you for your help!