       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: ClearAll@G

In: 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:  A = Table[G[n], {n, 4, 1000, 2}];

and plot A :

In:  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:  SMaxima[B_] := Module[{M}, {M = {}; AppendTo[M, B[]];
{Xg, Yg} = B[];
Do[If[Yg < B[[i + 1]]]],
{{Xg, Yg} = B[[i + 1]],
AppendTo[M, B[[i+ 1]]]}],
{i, 1, Length[B] - 1}]}; M]

In:  M = SMaxima[A]

and you get:

Out: 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:  plt1 = ListPlot[M, PlotJoined -> True, PlotRange -> All,
PlotStyle -> Hue[.6]]

Compare A and M :

In: Show[plt0, plt1]

Pick the Maximum point in M given by

{M[[Length[M]]][], M[[Length[M]]][]} :

In: Off[General::"spell1"]

In: Data = {{M[[Length[M]]][], M[[Length[M]]][]}}

Get:

Out: 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:  << Statistics`NonlinearFit`

In:  Clear[BF, alpha]

In:  BF[N_] = NonlinearFit[Data, alpha/(Exp[2/Log[N]] - 1),
{N}, {alpha}]

Out:  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: BoundSet = Table[{M[[i]][], BF[M[[i]][]]}, {i, 1, Length[M]}]

Get:

Out:

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: plt2 = ListPlot[BoundSet, PlotJoined -> True,
PlotRange -> All, PlotStyle -> Hue[.4]]

Plot A, M, and BoundSet together:

In:  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: Data2 = {M[], M[[Length[M]]]}

Get:

Out: Data2 = {{4, 1}, {990, 52}}

In: Clear[BF2, alpha, beta]

In: \$IterationLimit = Infinity

In: 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: BF2{N_] = 0.0187577*Exp[N^0.995686]

To complete our analysis we want to finish with:

In: BoundSet2 = Table[{M[[i]][], BF2[M[[i]][]]}, {i, 1, Length[M]}]

In: plt3 = ListPlot[BoundSet2, PlotJoined -> True,
PlotRange -> All, PlotStyle -> Hue[.4]]

and finally:

In: Show[plt0, plt1, plt3]

My questions are:

(1.) How can I improve the In 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?

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