Оценка максимального правдоподобия для обратного распределения Гаусса

Question

Я пытаюсь воспроизвести среднеквадратичную ошибку между фактическим и оценочным параметром 'tau' (for over a month :().Оценка 'tau', а именно 'tau_hat', получена путем оценки максимального правдоподобия (MLE), показанной ниже.

Совместная функция плотности вероятностиf(y|x,tau) задается как

, где u_i = x_i +T и T~IG(mu,lambda).И.Г .: Обратный гауссов.u - ожидаемое значение y.PDF f_T(t) задается как

Код, который я написал для этого веб-сайта , равен

    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

Результат, который я получаю:

Однако я должен получить тот, на который указывает красная стрелка

Какую ошибку я делаю в своем коде?Я не совсем уверен, как я рассчитал argmax.Для справки, научная статья, на которую я ссылаюсь, это здесь .

Answer 1

Я не могу запустить ваш код, так как для этого требуются некоторые наборы инструментов, которых у меня нет.Тем не менее, следующая строка:

tau_hat(ii)     =  max(Lkeh); 

даст вам значение максимума вероятности.Это не то, что вы действительно ищете, то есть tay_hat , при котором достигается ваше максимальное правдоподобие.

Вам нужна функция tay, которая отображает tay_hat в правдоподобие для данного значенияиз tay_hat.Предположим, это то, что вы делаете здесь, я не уверен, где зависимость от tay_hat.Предположим, что Lkeh - это то, что я только что описал, тогда

[maxLikelihoodValue, maxLikelihoodIndex] = max(Lkeh);

, используя оба выхода функции max, вы получите максимальное значение правдоподобия и, прежде всего, индекс, при котором этот максимум достигается.Если вы явно определили вектор tay, tay_hat будет дано просто как

tay_hat = tay (maxLikelihoodIndex);

, поэтому в основном это значение tay, для которого вы получаете максимальную вероятность, а несамо максимальное правдоподобие.

Чтобы дать вам игрушечный пример, предположим, что ваша функция правдоподобия L (x) = -x ^ 2 - 2 * x,

предположим, что она дискретна, так что

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

, тогда дискретная версия L будет

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

, тогда максимальное значение правдоподобия будет просто дано как

max(L_x);

, что составляет 0.9988 (на самом деле близко к точному значению)

, но вы получаете the value of x at which this maximum occurs.

Поэтому вы сначала извлекаете индекс в той последовательности, в которой вы получаете максимум, через:

[maximumLikelihood, maxLikIndex ] = max(L_x) ;

, а затем найдите оценку x по этому самому индексу, просто запрашивая значение x по этому индексу с:

x (maxLikIndex)

, что составляет около -1,0, как и ожидалось.В вашем примере вы хотите оценить наиболее вероятное значение tau_hat, которое (в структуре для часто используемых) определяется значением, которое максимизирует вашу функцию (которая не является максимумом самой функции).