fitting poisson

x = 0:3; % probable outcomes
freq = [36,40,22,2]; % frequencies of each outcome
n = 3; % number of trials
mean = x*freq'/sum(freq);
p=mean/n;
var=n*p*(1-p);

Px = binopdf(x,n,p); %P(X=x)
Exp = Px .* sum(freq); % expected frequencies

stem(x,freq);
hold;
plot(x,Exp,'--');
legend('given frequencies','expected frequencies');

figure();
plot(x,Px);
title("BINOMIAL PROBABILITY DISTRIBUTION");