matlab - 逆高斯分布的最大似然估计

Question

我正在尝试重现实际参数和估计参数 'tau' 之间的均方误差(一个多月 :()。估计的 ' tau'，即'tau_hat'是通过最大似然估计(MLE)得到的，如下所示。

联合概率密度函数f(y|x,tau)由

给出

其中 u_i = x_i +T 和 T~IG(mu,lambda)。 IG:逆高斯。 u 是 y 的期望值。f_T(t) 的 pdf 由

给出

我写的代码，基于这个website , 是

    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 的。供您引用，我指的科学论文是here .

Answer 1

我无法运行您的代码，因为它需要一些我没有的工具箱。也就是说，以下行:

tau_hat(ii)     =  max(Lkeh); 

将为您提供可能性的最大值。这不是您真正想要的，而是您达到最大可能性的 tay_hat 。

对于给定的 tay_hat 值，您需要一个 tay 函数将 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(实际上接近精确值)

但您所追求的是出现此最大值的 x 值。

因此，你首先提取你得到最大值的序列中的索引，通过:

[maximumLikelihood, maxLikIndex ] = max(L_x) ;

然后在该索引处找到 x 的估计值，只需使用以下命令请求该索引处的 x 值:

x (maxLikIndex)

正如预期的那样，大约是 -1.0。在您的示例中，您想要估计最有可能的 tau_hat(在常客框架中)由最大化函数的值(不是函数本身的最大值)给出。