function [nig_alpha,nig_beta,nig_mu,nig_delta]=nigem(sampoints,tol,maxiter)
  
  % NIGEM A Maximum Likelihood parameter estimation of the Normal-Inverse
  %  Gaussian using an EM type algorithm.
  %  
  %   [NIG_ALPHA, NIG_BETA, NIG_MU, NIG_DELTA] = NIGEM(SAMPOINTS, TOL, MAXITER) returns the scale 
  %   parameter NIG_ALPHA, the skewness parameter NIG_BETA, the location 
  %   parameter NIG_MU and the scale parameter NIG_DELTA, given the one
  %   dimensional set of points SAMPOINTS. The algorithm terminates if
  %   the number of iterations exceed MAXITER or 
  %   $max_{i=1..4} |\theta_i(k+1)-\theta_i(k)/theta_i(k)| < TOL$, where $i$
  %   is the iteration number and $\theta$ is the set of all parameters.
  %    
  %   SAMPOINTS = row vector, TOL = scalar, MAXITER = scalar
  %
  % 
  %   References: 
  %      [1] Karlis, D. (2002) An EM type algorithm for maximum
  %      likelihood estimation of the Normal-Inverse Gaussian distribution
  %
  % -----------------------------------------------------------------------
  % Saket Sathe
  %  email: saket@ee.iitb.ac.in
  % -----------------------------------------------------------------------
    
    
  %
  % General statistics of the input points.
  %
  samp_stddev = std(sampoints); % Standerd Deviation
  samp_var = var(sampoints);    % Variance, second moment
  samp_mean = mean(sampoints);  % Mean
  samp_kurt = kurtosis(sampoints); % Kurtosis, fourth moment
  samp_skew = skewness(sampoints); % Skewness, third moment

  %
  % Initial estimates of alpha, beta, mu, delta and gamma.
  %
  nig_gamma1 = samp_skew/((samp_var)^(3/2)); % gamma_1
  nig_gamma2 = samp_kurt/((samp_var^2)-3); % gamma_2
  nig_gamma = 3/(samp_stddev*sqrt(3*nig_gamma2-5*(nig_gamma1^2))) % gamma
  nig_beta = (nig_gamma1*samp_stddev*(nig_gamma^2))/3 % beta
  nig_delta = (samp_var*(nig_gamma^3))/(nig_beta^2+nig_gamma^2) % delta
  nig_mu = samp_mean-((nig_beta*nig_delta)/(nig_gamma)) % mu
  nig_alpha = sqrt(nig_gamma^+nig_beta^2) % alpha


  for iter=1:maxiter
    disp(sprintf('Iteration: %f',iter));

    %
    % E-step
    %
    
    s_i = zeros(size(sampoints));
    w_i = zeros(size(sampoints));

    for i=1:size(sampoints,2)
      temp_phi = (1 + ((sampoints(i)-nig_mu)/nig_delta)^2)^0.5;
            
      if besselk(1,(nig_delta*nig_alpha*temp_phi))~=0 && nig_delta*nig_alpha*temp_phi>0
        s_i(1,i) = (nig_delta*temp_phi*besselk(0,(nig_delta*nig_alpha*temp_phi)))/(nig_alpha*besselk(1,(nig_delta*nig_alpha*temp_phi)));    
      else
        disp('zero')
        nig_delta*nig_alpha*temp_phi
      end
  
      if besselk(-1,(nig_delta*nig_alpha*temp_phi))~=0 && (nig_delta* ...
                                                      nig_alpha*temp_phi) >0
        w_i(1,i) = (nig_alpha*besselk(-2,(nig_delta*nig_alpha*temp_phi)))/(nig_delta*temp_phi*besselk(-1,(nig_delta*nig_alpha*temp_phi)));
      else
        disp('zero')
        nig_delta*nig_alpha*temp_phi
      end
      
    end
    % nig_gamma
    old_params = [nig_beta nig_delta nig_mu nig_alpha];

    % M-step

    M_hat = sum(s_i,2)/size(s_i,2);
    %lamda_cap_hat = size(s_i,2)*(sum(((w_i-(M_hat^-1)),2))^-1;

    tempvar = w_i-(1/M_hat);

    lamda_cap_hat = size(s_i,2)/(sum(tempvar,2));

    nig_delta = lamda_cap_hat^0.5;

    nig_gamma = nig_delta/M_hat;

    nig_beta = (dot(sampoints,w_i) - samp_mean*sum(w_i,2))/(size(s_i,2)- ...
                                                  M_hat*sum(w_i,2));

    nig_mu = samp_mean - nig_beta*M_hat;

    nig_alpha = (nig_gamma^2+nig_beta^2 )^0.5;

    % nig_gamma
    new_params=[nig_beta nig_delta nig_mu nig_alpha];

    paramdiff = new_params-old_params;
    perdiff = paramdiff./old_params;

    if max(abs(perdiff))<tol
      nig_alpha
      nig_beta
      nig_mu
      nig_delta
      break
    end

end