Mathematica program to obtain a bounding function for a set of data points.
- To: mathgroup at smc.vnet.net
- Subject: [mg43248] Mathematica program to obtain a bounding function for a set of data points.
- From: gilmar.rodriguez at nwfwmd.state.fl.us (Gilmar Rodríguez Pierluissi)
- Date: Sat, 23 Aug 2003 08:08:23 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
The following is to help clarify my question about how to build a Mathematica program to obtain a bounding function for a set of data points (please refer to: http://forums.wolfram.com/mathgroup/archive/2003/Jul/msg00496.html). I haven't received a response so far, and I figure that the reason for this is that I need to elaborate my question better. (** Preliminary (but, necessary)information :**) Remark: Q_i below means "Q subscript i", etc. Let N be an integer greater or equal to 4. Assume the existance of primes P_i in the interval [2, N/2] in X, and primes Q_i = (N - P_i) in the interval [N/2, N - 2] in Y, with P_i <= Q_i such that N = P_i + Q_i. Then the set of points(P_i,Q_i) is called a Prime Partition of N. The Minimal Goldbach Decomposition Point ("MGDP" for short) is a point (P,Q) among the points (P_i,Q_i) such that the perpendicular distance between the points (P_i,Q_i) (located on the line Y = -X + N) and the point (N/2,N/2)(located on the line Y=X) is the shortest distance, i.e. perpdist[(P,Q), (N/2,N/2)]=Min[Sqrt[(P_i - N/2)^2+(Q_i- N/2)^2]] The following web pages help to clarify the concepts considered so far : http : www.gilmarlily.net goldbach gcstated.htm http : www.gilmarlily.net goldbach geoview.htm http : www.gilmarlily.net goldbach mgdplot1.htm (** End of preliminary stuff. **) First step; the following module gives the Minimal Goldbach Decomposition for an integer n >= 4. In[1]: MGDP[n_] := Module[{p, q}, {m = n 2; If[(m \[Element] Primes), {(p = m), (q = m)}, {k = PrimePi[m]; Do[If[(n - Prime[i]) \[Element] Primes, {hit = i, Break[]}], {i, k, 1, -1}], p = Prime[hit], q = (n - p)}]}; {p, q}] For example; In[2]: MGDP[100] Out[2]: {47, 53} Second step; we build a Module to calculate the perpendicular distance between the MGDP and the point (n/2, n/2). For each n, the Module returns a value {n, perpendicular distance}. In[3]: perpdist[n_] := Module[{p, q, d}, {m = n 2; dist[{x1_, y1_}, {x2_, y2_}] := Sqrt[(x2 - x1)^2 + (y2 - y1)^2]; If[(m \[Element] Primes), {(p = m), (q = m), (d = 0)}, {k = PrimePi[m]; Do[If[(n - Prime[i]) \[Element] Primes, {advpt = i, Break[]}], {i, k, 1, -1}], p = Prime[advpt], q = (n - p), d = dist[{p, q}, {m, m}]}]}; {n, d}] For example; In[4]: perpdist[100] Out[4]: {100, 3Sqrt[2]} Third step; for each n between 4 and 1 Million, we calculate the perpendicular distance between the MGDP of n, and the point (n/2,n/2). (The calculation of the set A might take a little time): In[5]: A = Table[perpdist[n], {n, 4, 10^6, 2}]; Fourth step; we seek to find a Bounding Function for the set A above. To start, we first build a module called SMaxima[]. This module compares each pair {(a_x,a_y), (a_(x+1),a_(y+1))} in A, and appends (a_(x+1),a_(y+1)) to a set M, if a_y < a_(y+1). In[6]: SMaxima[B_] := Module[{M}, {M = {}; AppendTo[M, B[[1]]]; {Xg, Yg} = B[[1]]; Do[If[Yg < B[[i + 1]][[2]], {{Xg, Yg} = B[[i + 1]], AppendTo[M, B[[i + 1]]]}], {i, 1, Length[B] - 1}]}; M] Applying SMaxima[] to A gives us : In[7]: M = SMaxima[A] Out[7]: {{4, 0}, {8, Sqrt[2]}, {16, 3Sqrt[2]}, {44, 9Sqrt[2]}, {92, 15Sqrt[2]}, {242, 18Sqrt[2]}, {256, 21Sqrt[2]}, {272, 27Sqrt[2]}, {292, 33Sqrt[2]}, {476, 39Sqrt[2]}, {530, 42Sqrt[2]}, {572, 45Sqrt[2]}, {682, 48Sqrt[2]}, {688, 75Sqrt[2]}, {1052, 87Sqrt[2]}, {1808, 93Sqrt[2]}, {2228, 117Sqrt[2]}, {3382, 120Sqrt[2]}, {3472, 135Sqrt[2]}, {3502, 138Sqrt[2]}, {3562, 168Sqrt[2]}, {4952, 183Sqrt[2]}, {6194, 210Sqrt[2]}, {7102, 228Sqrt[2]}, {10262, 300Sqrt[2]}, {17008, 333Sqrt[2]}, {20684, 369Sqrt[2]}, {37052, 393Sqrt[2]}, {45128, 453Sqrt[2]}, {49552, 525Sqrt[2]}, {80144, 621Sqrt[2]}, {137414, 720Sqrt[2]}, {251806, 810Sqrt[2]}, {349826, 846Sqrt[2]}, {362534, 1086Sqrt[2]}, {742856, 1281Sqrt[2]}} Plot M using the command: In[8]: plt1 = ListPlot[M, PlotJoined -> True, PlotRange -> All, PlotStyle -> Hue[.6]] Out[8]: You get a plot. Fifth step; the set M becomes instrumental in finding a Bounding Function for A : In[9]: << Statistics`NonlinearFit` We want our bounding function to pass through the point in M with the highest y - coordinate, so we do : In[10]: Off[General::"spell1"] In[11]: data = {{M[[Length[M]]][[1]], M[[Length[M]]][[2]]}} Out[11]: {{742856, 1281Sqrt[2]}} In[12]: BF[N_] = NonlinearFit[data, Alpha/(Exp[2/ Log[N]] - 1), {N}, {Alpha}] Out[12]: (288.8656741801377)/(-1 + Exp[2/Log[N]) We now build a set called "Bounding Set" (i.e. "BoundSet" for short) where the points have the same x - coordinates as the points in M, but the y - coordinates are calculated using our Bounding Function above : In[13]: BoundSet = Table[{M[[i]][[1]], BF[M[[i]][[1]]]}, {i, 1, Length[M]}] Out[13]: {{4, 89.3744}, {8, 178.71}, {16, 273.236}, {44, 414.792}, {92, 519.275}, {242,657.102}, {256, 665.136}, {272, 673.799}, {292, 683.94}, {476, 753.851}, {530, 769.239}, {572, 780.162}, {682, 805.363}, {688, 806.618}, {1052, 867.505}, {1808, 945.222}, {2228, 975.218}, {3382, 1035.18}, {3472, 1038.95}, {3502, 1040.19}, {3562, 1042.63}, {4952, 1089.99}, {6194, 1122.17}, {7102, 1141.84}, {10262, 1194.79}, {17008, 1267.49}, {20684, 1295.65}, {37052, 1379.59}, {45128, 1407.98}, {49552, 1421.45}, {80144, 1490.7}, {137414, 1568.38}, {251806, 1655.66}, {349826, 1703.05}, {362534, 1708.19}, {742856, 1811.61}} We now plot the points in BoundSet : In[14]: plt2 = ListPlot[BoundSet, PlotJoined -> True, PlotRange -> All, PlotStyle -> Hue[.4]] Out[14]: You get a second plot. Finally, we plot the sets M and BoundSet together : In[15]: Show[plt1, plt2] Out[15]: You get a third plot. The point that I'm really trying to make via this example is : it would be nice if Wolfram Research could combine the "SMaxima[] technique" above, together with the NonlinearFit program to come up with a new program that calculates the Bounding Function for an arbitrary set of points similar to the set A above. Thank you for your attention and your help!