4
$\begingroup$

"Assume independent random variables $Y_i$~$Poisson(λx_i)$. Supposing that $x_i$ are given, fixed constants, obtain the least squares estimator of $λ$ and compute its variance."

This kind of a problem is the first one faced by me. The first things that are coming to my mind are regression equation of the form $log(Y)=a+bx$ and then trying the least square algorithm by minimizing $∑(log(y)-a-bx)^2$. From this I will be getting the least squares estimates of $a,b$ but how do I get $λ$?

As I said, this is the first time I am facing such a question. I know point estimation in terms of MLE and moment-equality. I need help and resources (i.e. more knowledge) about solving this kind of a question. Can you kindly supply me with the answer and the explanation along with maybe links to some pdf's? Google is not helping me. Thanks in advance.

(EDIT:) While I was waiting for any response to this question, I looked up deeper and found that I did not know one information. That is, the least square estimator for a parameter $θ$ is obtained by minimising $∑(Y_i-E(Y_i|x_i,θ))^2$. So I tried to work on this using my new knowledge.

I get the LSE of $λ$ as $\dfrac{∑x_iY_i}{∑(x_i)^2}$. Please tell me if this is correct.

Then the variance of this LSE is $\dfrac{λ∑(x_i)^3}{(∑(x_i)^2)^2}$.

Is this right? Thanks for your time.

$\endgroup$
  • 2
    $\begingroup$ If $Y_i \sim \mbox{Poisson}(\lambda x_i)$, what is $E(Y_i | x_i)$? Could you minimize $\sum \left( Y_i - E(Y_i | x_i)\right)^2$? As far as taking the log goes, what happens if $Y_i =0$, which occurs with positive probability under a Poisson distribution? $\endgroup$ – Andrew M Oct 31 '14 at 7:22
  • $\begingroup$ Dear Andrew, thank you for this information. I did not initially know that the true method is minimizing $∑(Y_i−E(Y_i|x_i))^2$. Fortunately, I searched it out myself on Google shortly before you commented. Thank you. Would you kindly verify if my answers are correct? I have edited my answer. $\endgroup$ – Landon Carter Oct 31 '14 at 7:27
3
$\begingroup$

Log Likelihood Solution

The Log Likelihood Function is given by:

$$ \begin{align*} \hat{\lambda} & = \arg \max_{\lambda} p \left( \boldsymbol{y} \mid \boldsymbol{x}, \lambda \right) \\ & = \arg \max_{\lambda} \log \prod_{i = 1}^{n} \frac{ \left( \lambda {x}_{i} \right)^{ {y}_{i} } }{ {}y_{i}! } {e}^{-\lambda {x}_{i}} \\ & = \arg \max_{\lambda} \sum_{i = 1}^{n} \log \left( \frac{ \left( \lambda {x}_{i} \right)^{ {y}_{i} } }{ {}y_{i}! } {e}^{-\lambda {x}_{i}} \right) \\ & = \arg \max_{\lambda} \log \lambda \sum_{i = 1}^{n} {y}_{i} + \sum_{i = 1}^{n} \log {x}_{i} - \lambda \sum_{i = 1}^{n} {x}_{i} -\sum_{i = 1}^{n} \log {y}_{i} ! \\ & \Rightarrow \hat{\lambda} = \frac{ \sum_{i}^{n} {y}_{i} }{ \sum_{i}^{n} {x}_{i} } \end{align*} $$

This indeed coincide with the classic case where $ {x}_{i} = 1 $ and then the MLE is the empirical average.

Least Squares Solution

The Least Squares solution is given by:

$$ \begin{align*} \hat{\lambda} & = \arg \min_{\lambda} \sum_{i}^{n} \left( \mathbb{E} \left[ {y}_{i} \mid {x}_{i}, \lambda \right] - {y}_{i} \right)^{2} \\ & = \arg \min_{\lambda} \sum_{i}^{n} \left( \lambda {x}_{i} - {y}_{i} \right)^{2} \\ & \Rightarrow \hat{\lambda} = \frac{ \sum_{i}^{n} {x}_{i} {y}_{i} }{ \sum_{i = 1}^{n} {x}_{i}^{2} } \end{align*} $$

Small simulation in MATLAB:

% Mathematics Q122153
% https://stats.stackexchange.com/questions/122153
% Least Squares Estimation of Poisson Parameter
% References:
%   1.  Poisson Distribution Wikipedia - https://en.wikipedia.org/wiki/Poisson_distribution.
%   2.  Poisson Distribution (MATLAB)- https://www.mathworks.com/help/stats/poisson-distribution.html.
%   3.  Poisson Random Numbers (MATLAB)- https://www.mathworks.com/help/stats/poissrnd.html.
% Remarks:
%   1.  sa
% TODO:
%   1.  ds
% Release Notes
% - 1.0.000     10/09/2017
%   *   First release.


%% General Parameters

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = OFF;


%% Simulation Parameters

paramLambda = 0.75;
numSamples  = 1000;


%% Algorithm Parameters

gridStartVal    = 0.01;
gridEndVal      = 2;
numGridSamples  = 2000;


%% Generate Data

vX              = randi([1, 10], [numSamples, 1]); %<! Known
vParamLambda    = paramLambda * vX;
vDataSamples    = poissrnd(vParamLambda, [numSamples, 1]);

vLambdaGrid     = linspace(gridStartVal, gridEndVal, numGridSamples);


%% Maximum Likelihood Estimator of Lambda - Brute Force

vLogLikelihood  = zeros([numGridSamples, 1]);

for ii = 1:numGridSamples
    currLambda = vLambdaGrid(ii);
    vLogLikelihood(ii) = log(currLambda) * sum(vDataSamples) + sum(log(vX)) - currLambda * sum(vX) - sum(log(factorial(vDataSamples)));
end


%% Maximum Likelihood Estimator of Lambda - Closed Form

paramLambdaMle = sum(vDataSamples) / sum(vX);


%% Least Squares Estimator of Lambda - Brute Force

vLsLikelihood  = zeros([numGridSamples, 1]);

for ii = 1:numGridSamples
    currLambda          = vLambdaGrid(ii);
    vLsLikelihood(ii)   = sum( (vDataSamples - (currLambda * vX)) .^ 2 );
end


%% Least Squares Estimator of Lambda - Closed Form

paramLambdaLs = sum(vDataSamples .* vX) / sum(vX .^ 2);


%% Analysis

hFigure     = figure('Position', figPosLarge);
hAxes       = subplot(2, 1, 1);
set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vLambdaGrid, vLogLikelihood);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries = plot([paramLambda, paramLambda], get(hAxes, 'YLim'));
set(hLineSeries, 'LineStyle', ':');
hLineSeries = plot([paramLambdaMle, paramLambdaMle], get(hAxes, 'YLim'));
set(hLineSeries, 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['Log Likelihood Objective Function Value vs. \lambda Hypothesis']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', '\lambda', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Objective Function Value', ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['Log Likelihood Objective Function'], ['Log Likelihood Maximum'], ['Ground Truth']});
% set(hAxes, 'LooseInset', [0.07, 0.07, 0.07, 0.07]);

hAxes       = subplot(2, 1, 2);
set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vLambdaGrid, log(vLsLikelihood));
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries = plot([paramLambda, paramLambda], get(hAxes, 'YLim'));
set(hLineSeries, 'LineStyle', ':');
hLineSeries = plot([paramLambdaLs, paramLambdaLs], get(hAxes, 'YLim'));
set(hLineSeries, 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['Log Least Squares Objective Function Value vs. \lambda Hypothesis']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', '\lambda', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Objective Function Value', ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['Log Least Squares Objective Function'], ['Log Least Squares Minimum'], ['Ground Truth']});
% set(hAxes, 'LooseInset', [0.07, 0.07, 0.07, 0.07]);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end


%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

enter image description here

The full code is available on my StackExchange Cross Validated Q122153 GitHub Repository.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.