function [a,b] = gam_wfit( x, w, nsamples )
  
  if nargin < 2
    error( [ mfilename ' needs 2 input arguments.' ] );
  end
  
 if ~all( size( x ) == size( w ) )
   error( [ mfilename ': "x" and "w" must have the exact same size.' ] );
 end

 if any( isnan( x(:) ) | isinf( x(:) ) )
   error( [ mfilename ': "x" must not contain any Inf or NaN value.' ] );
 end

 if any( isnan( w(:) ) | isinf( w(:) ) )
   error( [ mfilename ': "w" must not contain any Inf or NaN value.' ] );
 end

 if ~exist( 'nsamples', 'var' )
   nsamples = 1000;
 end
 
 % Now we can drop size information
 x = x(:);
 w    = w(:);
 
 % Order x by increasing value
 [ x, ind ] = sort( x );
 w = w( ind );
 
 % Compute the CDF
 cumw = cumsum( w );
 cumw = cumw / cumw( end );
 
 [ cumw, ind ] = unique( cumw );
 cumx          = x( ind );
 
 % Invert it: find 1000 representative samples
 
 samples = interp1( cumw, cumx, (0.5:nsamples)/nsamples, 'linear', 'extrap' );
 
 % Fit a gamma on this
 phat = gamfit( samples( find( samples > 0 ) ) );
 
 % Return the parameters
 a = phat( 1 );
 b = phat( 2 );