function [alpha,r_squared] = PLmaxlikelihood(x,xmin)

if (nargin < 2)
    xmin = 1;
end
if ((min(x) <= 0)&(nargin < 2))
      disp('warning: x cannot be 0 or negative, will fit only for x > 0');
      x = x(x>0);
end

if (xmin<= 5)
      disp('warning: fit may not be accurate for xmin < 5 for discrete data');
end

xx = x(x >= xmin);
mysum=0;
for i=1:length(xx)
    mysum = mysum+log(xx(i)/xmin);
end

alpha = 1+length(xx)/mysum;

mysum = 0

xfit = min(x):max(x);
for i=min(x):max(x)
    mysum = mysum + 1/i^alpha;
end

[a,b] = cumulativecounts(x);
b = b/b(1);
[aa,bb] = binvector(x);
bb = bb/sum(bb);
fitline = min(x):max(x);
fitline = 1/mysum*1./xfit.^alpha;
for i=1:length(fitline)
    fitlinecumulative(i) = sum(fitline(i:length(fitline)));
end
%sum(fitline)

size(bb)
size(fitline)
figure(1);
p = loglog(xfit(1:length(a)),fitlinecumulative(1:length(a)),'b-',a,b,'rs:')
legend('fit','data')
axis([min(a) max(a) min(b) max(b)])
xlabel('x')
ylabel('cumulative distribution')

rsquared = 1-sum((bb - fitline).^2)/sum((bb - mean(bb)).^2)
alphastr = sprintf('alpha = %.4f fit, R^2 = %.6f',alpha,rsquared);

figure(1);
p = loglog(xfit(1:length(a)),fitline(1:length(a)),'b-',aa,bb,'rs')

legend(alphastr,'data')
set(p,'MarkerSize',2);
axis([min(a) max(a) min(b) max(b)])
xlabel('x')
ylabel('histogram')

figure(2);
p = loglog(xfit,fitlinecumulative,'b-',a,b,'rs')
legend(alphastr,'data')
set(p,'MarkerSize',2);
axis([min(a) max(a) min(b) max(b)])
xlabel('x')
ylabel('cumulative histogram')

diffvec = bb-fitline(1:length(bb));
xx = 1:length(b);

%var(bb)
%semilogx(xx,diffvec,'b-',xx,bb,'r-')