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]?