more off
				# random data
clear all;
N=2000; # if this is changed then may need to alter intendedk
M=5;   # true dimension of space
intendedk = [ 2 3 4 6 8 16 32 64 128 256 ] ;
[junk,sR] = size(R) ;

x = rand( N, M ) ;
sitot= zeros(size(R)) ;
sicount = zeros(size(R)) ;
numerator = zeros( size(R) ) ;   ##
denominator = zeros( size(R) ) ; ## for the correct estimator

for n=1:N
  d = x - ones(N,1)*x(n,:) ; # the matrix of all distances to the nth point

  d = d + (d>0.5)*(-1) ;
  d = d + (d<-0.5)*(1) ;

  T = sum ( (d.^2)'  )'  ; ## sum over horizontal dimension
  S = sort(T) ;

  lT = log(S) * 0.5 ;
  lT(n) = 0;

  for r=1:sR
    k = intendedk(r) ;
    slT = sum( lT(2:k) ) ;
    mlT = slT / (k-1) ;
    lR  = lT(k+1) ;
    imhat(n,r) = lR-mlT ;
    mhat(n,r) = 1.0/(lR-mlT) ;

    numerator(r) = numerator(r) + ( (k-1) * lR - slT ) ;
    denominator(r) = denominator(r) + (k-1) ;

  endfor
end

theirestimator = mean( mhat ) ;

imhatmean = mean( imhat ) ;

anotherm = 1.0./imhatmean ;
				# another and my should be the same (above and below)
myestimator = numerator ./ denominator
mym = 1.0./myestimator

## mhat = 1./imhat ;

doplots
