CODEWITHSUNDEEP
matlab
CasulOrnstein
CasulOrnstein

Reputation: 23

Fitting a upper bound to a set of points

I have the following code about the Gauss circle problem, which finds the error between the number of lattice points of circles, of varying radius R, and their area. It then plots a graph of R against the error.

  %LatticeError plots graph of error between the lattice points in a circle of radius R and pi R^2, against R. 
R = 0:0.5:500;
  %Calculates the area of each circle
Area = pi .* (R.^2);
  %Creates an array of the number of points contained in each circle
LatticePts = arrayfun(@lattice,R); 
Error = abs(LatticePts - Area);
plot(R,Error)

What I want to be able to do is be able to use Matlab to approximate a function, f(R) = a*(R^b), which is an upper bound to the error over the range of R. But I'm not sure how to do so, any ideas?

Upvotes: 2

Views: 1526

Answers (1)

user3717023
user3717023

Reputation:

Putting a curve of the form y=ax^b over a given set of points (xi,yi) amounts to putting a line over the set (log(xi), log(yi)), since a power function becomes linear on the log-log scale.

Applying convhull to the set of (log(xi), log(yi)) gives a convex hull of these points. Any edge of this convex hull gives a line that fits tightly against the data set. Strictly speaking, there is no "best" line among these, but it's reasonable to pick one that is lowest in the middle of your range of values (that is around R=250). This is what I did below. For completeness, I included the computation of the number of lattice points at the beginning.

Computation of Error

Radius = [];
Error = [];
for R = 0.5:0.5:500;
    Area = pi*R^2;
    n = floor(R);
    x = -n:n;
    y = -n:n;
    [X,Y] = meshgrid(x,y);
    LatticePts = nnz(X.^2+Y.^2<=R^2);
    Radius = [Radius, R];
    Error = [Error, max(1,abs(LatticePts - Area))];
end

The reason to take maximum with 1 is to avoid any problems with logarithms of zero or near-zero numbers.

Convex hull on log-log scale

ind = convhull(log(Radius), log(Error));

This gives the indices of extreme points, listed counterclockwise. This is actually the opposite order from the one I want: it begins with lower bounds, which are uninteresting. The ones on top are at the end of the list, so I search through it starting at the end:

Side of convex hull lying over R=250:

i = numel(ind);
while (Radius(ind(i))<250)
    i = i-1;
end

Computation of a,b and output

R1 = Radius(ind(i));
E1 = Error(ind(i));
R2 = Radius(ind(i+1));
E2 = Error(ind(i+1));
b = log(E1/E2)/log(R1/R2);
a = E1/R1^b;
plot(Radius, Error, '.');
hold on
plot(Radius, a*Radius.^b , 'r-');
hold off
fprintf('a = %.3f, b = %.3f\n', a,b);

I got a = 3.219, b = 0.617. Here is the plot:

gauss

Upvotes: 2

Related Questions