Maximum Likelihood estimation for Inverse Gaussian distribution

Question

I am trying to reproduce the mean squared error between the actual and estimated parameter 'tau' (for over a month :(). The estimated 'tau', namely 'tau_hat' is obtained through the maximum likelihood estimation (MLE) , shown below.

The joint probability density function f(y|x,tau) is given by

where u_i = x_i +T and T~IG(mu,lambda). IG: Inverse Gaussian. u is the expected value of y. The pdf of f_T(t) is given by

The code I have written, based on this website, is

    clear
    lambda  =   8.1955;
    mu      =   10;
    N       =   128; % max number of molecules
    x       =   zeros(N,1); % transmission time of the molecules from the Tx; for K = 1
    tau     =   .5; % arbitrary initital tau
    simN    =   1000 ; % # runs per N 
    no_molecules_per_simN   =  [4, 8, 32, 64, N];
    tau_hat   =   zeros(size(no_molecules_per_simN));

    for ii=1: length(no_molecules_per_simN)

        Lkeh  = zeros(1,length(no_molecules_per_simN(ii)));  % inititalize likelihood array

        for jj=1: simN
            T               =  random('InverseGaussian', mu,lambda, [no_molecules_per_simN(ii),1]); % random delay
            y_prime         =  x(1:no_molecules_per_simN(ii)) + T + tau; % arrival time of the molecules seen by the Rx
            y_prime_sort    =  sort(y_prime); % to arrange them in the ascending order of arrival
            u               =  y_prime_sort;  % assign to u variable
            t               =  u - x(1:no_molecules_per_simN(ii)) - tau;
            for kk = 1: length(u)
                % applying the likelihood function to eq. 3 and ignoring the constant terms
                 %linear likelihood
%             Lkeh(jj,kk)    =  prod(t(kk).^-1.5).*exp(-sum((t(kk) - mean(t)).^2./t(kk)).*(lambda./(2.*mean(t).^2 )));

% [UPDATE to the code]
            % log likelihood
            Lkeh(jj,kk)    =   -1.5*sum(t(kk))-(lambda./(2.*mu.^2 )).*sum((t(kk) - mu).^2./t(kk));

            end

        end
        Lkeh_mean       =  mean(Lkeh,1); % averging the values
    % [UPDATE to the code]
        [maxL,index]    =  max(Lkeh_mean);
        t_hat(ii)       =   T(index) ; % this will give the likelihood value of the propagation delay
        tau_hat(ii)     =   mean(u -  x(1:no_molecules_per_simN(ii)) - t_hat(ii)); % reverse substitution

    end

    MSE = zeros(size(tau_hat)); % initializing the array for MSE

    for ii=1:length(tau_hat)
        MSE(ii) = immse(tau,tau_hat(ii)); % mean squared error
    end

    figure
    loglog(no_molecules_per_simN,MSE,'-o')
    xlabel('n_{1}(quantity of molecules)')
    ylabel('MSE(sec^{2})')
    grid on

The result I obtain is

However, I should be obtaining the one pointed to by the red arrow

What is the mistake I am making in my code? I am not quite sure of how I calculated the argmax. For your reference, the scientific paper I am referring to is here.

Answer 1

I cannot run your code, since it requires some toolboxes, which I do not have. That said, the following line :

tau_hat(ii)     =  max(Lkeh); 

will give you the value of the max of the likelihood. That is not what you are really after, which is the tay_hat at which your maximum likelihood is attained.

You need a function of tay which maps tay_hat into the likelihood, for a given value of tay_hat. Suppose that is what you do here, I am not sure where the dependency on tay_hat is. Let assume Lkeh is what I have just described then

[maxLikelihoodValue, maxLikelihoodIndex] = max(Lkeh);

using both outputs of the max function, you will have the maximum likelihood value and, above all, the index at which that max occurs. If you have explicitly defined the tay vector the tay_hat will be given simply by

tay_hat = tay (maxLikelihoodIndex) ;

so basically is the value of tay for which you get the maximum likelihood, not the maximum likelihood itself.

To give you a toy example, suppose your likelihood function is L(x) = -x^2 - 2*x,

suppose it is discretised so that

 x = linspace(-2,2,30);

then the discrete version of L will be

L_x = -x.^2 -2*x;

then the maximum likelihood value will be simply given by

max(L_x);

which happens to be 0.9988 (in fact close to the exact value)

but what you are after is the value of x at which this maximum occurs.

Therefore, you first extract the index in the sequence at which you get the maximum, through:

[maximumLikelihood, maxLikIndex ] = max(L_x) ;

and then find the estimate of x at that very index simply requesting the value of x at that index with :

x (maxLikIndex)

which is about -1.0, as expected. In your example, you want to estimate your most likely tau_hat which (in a frequentist framework) is given by the value which maximises your function (which is not the maximum of the function itself).