% demo on score function and Fisher information
% this builds up example 3.4.2 in Bickel & Doksum, Mathematical Statistics vol. 1 (2nd ed)
% Jerry Zhu

% pdf of Poisson P(k; theta) = theta^k exp(-theta) / k!
%theta=1 % true parameter
theta=3 % true parameter

hold on
for trial=1:10000,

	% a sample X1...Xn ~ Poisson(theta)
	n = 10;
	for i=1:n
	  X(i) = poissrnd(theta);
	end

	% loglikelihood function Lambda(theta) = log P(X;theta) = sum_{i=1}^n [X(i)log(theta) - theta - log(X(i)!)]
	th=0.01:0.01:6;
	Lambda = sum(X)*log(th) - n*th - sum(log(factorial(X)));
	plot(th, Lambda, 'b-')

	% the score function s(theta,X)=
	% partial log(P(X;theta)) / partial theta  = sum_{i=1}^n [ X(i)/theta - 1] = sum(X)/theta - n
	s(trial) = sum(X)/theta - n;

	% the MLE on the sample
	MLE(trial) = sum(X)/n;


end

% we expect E[s]=0
figure
hist(s,100)
mean(s)

% Fisher info = var(s) = var(sum(X)/theta - n) 
% note sum(X) is itself a Poisson with rate n*theta, so
% var(sum(X)/theta - n) = var(sum(X)/theta) = var(sum(X)) / theta^2 = n*theta / theta^2 = n/theta
I = n/theta
% Monte Carlo
var(s)

mean(MLE)
var(MLE)
