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

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)的概率密度函数如下：

我根据这个网站编写的代码是：

    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。供您参考，我所参考的科学论文在这里。

- nashynash

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（实际上接近于精确值）。

但你需要的是 使得该最大值发生的x的值。

因此，首先通过以下方式提取序列中达到最大值的索引：

[maximumLikelihood, maxLikIndex ] = max(L_x) ;

然后，在该索引处找到x的估计值，只需使用以下方法请求该索引处的x值：

x (maxLikIndex)

这个值约为-1.0，符合预期。 在你的例子中，你想要估计出最可能的tau_hat值，在频率学框架中，这个值由使函数取最大值（而不是函数本身的最大值）的数值给出。

- pacta_sunt_servanda

---

我已经尝试过了，但是没有得到想要的结果。 - nashynash

我不是你领域的专家，但你需要的是数据来计算可能性；你是通过抽样生成合成数据吗？在我看来似乎是这样。如果我已经理解了你的代码，我会添加一些必要的点来干扰你的代码。 - pacta_sunt_servanda

为了数值稳定性，使用对数似然和似然本身是方便的（因此您的prod将需要被替换为sum）。 - pacta_sunt_servanda

2. u是y_prime_sort，已经排序过了。我不确定为什么需要排序，但是假设你想要使用排序后的向量进行操作，那么你需要对所有向量应用相同的排序，否则这会影响你的结论。同样，你可能想要使用逻辑索引；[y_prime_sort，sortingOrderIdx] = sort(y_prime); 然后请将sortingOrderIdx应用到你在计算中使用的其他向量上。

- pacta_sunt_servanda

1. 我也使用了对数似然函数。[代码也将进行更新]
2. 'u' 被排序以按照它们到达的时间顺序排列分子。

- nashynash