Poisson distribution

The Poisson distribution is a discrete probability expressing the probability of a given number of events occurring in a fixed interval where the events occur with a known average independent rate, λ .

p(k) = λ ^k exp(-λ ) /factorial(k)

set.seed(125)
n<-1000
lambda<-5
poik<-rpois(n, lambda)
barplot(table(poik)/n,col=3,main="Poisson distribution n=1000")

poi20<-rpois(20, lambda)
barplot(table(poi20)/10,col=4,main="Categolical distribution n=20")

Its probability for the occurence in a fixed interval is λ .
We can make function named as PrPoi(x, lambda), where x is a vector of numbers of which the events occur in that interval.
This function will back the value of
Π {exp(-λ )(λ ^x[i])/x[i]! }

PrPoi<-function(x,lambda){
  p<-1
  for(ki in 0:max(x)){
    p<-p*exp(-lambda)*lambda^ki/factorial(ki)
  }
  p
}
#
PrPoi(x=poik,lambda=lambda)
## [1] 1.533356e-18
PrPoi(x=poi20,lambda=lambda)
## [1] 2.987485e-12

Likelifood distribution & Maximum likelifood (ML)

Distribution of λ in the world may be estimated and visualized as a graph from the data, x. Maximum likelihood of world’s λ may be gained also from this data distribution.
Lambda parameter as the axis of the world are 1 =< λ < = n. Where n =< max(x)
Along this axies, probabilities of the likelifood will be estimated.

Likelifood.Poi<-function(x){ 
  Lk<-NULL
  for( lam in 0:max(x)){
     p<-PrPoi(x=x,lambda=lam)
     Lk<-rbind(Lk,c(lam,p))
  }
  Lk[,2]<-Lk[,2]/sum(Lk[,2])
  Lk<-data.frame(Lk)
  names(Lk)<-c("lambda","Likelihood")
  Lk
}
#####
#n=1000
LkPoi<-Likelifood.Poi(x=poik)
plot(LkPoi,col=3,type="l",main="Likelihood")
#ML
LkPoi[LkPoi$Likelihood==max(LkPoi$Likelihood),]
##   lambda Likelihood
## 7      6  0.5867565
#n=10
LkPoi<-Likelifood.Poi(x=poi20)
lines(LkPoi,col=4,type="l")

#ML
LkPoi[LkPoi$Likelihood==max(LkPoi$Likelihood),]
##   lambda Likelihood
## 6      5  0.4579224

Gamma Distribution

The gamma distribution is a two-parameter family of continuous probability distributions. In R, these are shape and rate.
This is faverable in Bayesian statistics, where the gamma distribution is used as a conjugate prior distribution for various types of inverse scale parameters, such as the λ of an exponential distribution or a Poisson distribution.

p(y) = rate^shape/gamma(shape) y^(shape-1)exp(-rate y) = y^(shape-1)exp(-y)/gamma(shape) (iff. rate = 1)

m<-max(c(poi20,poik))
prior<-dgamma(x=0:m,shape=5,rate=1)
plot(x=0:m,y=prior,col=2,type="l", main="Gamma distribution",xlab="lambda")

#Max
which(prior==max(prior))-1
## [1] 4

Conjugate Distributions

As you see, both have the same form
x^a exp(-x)
So
PrPoi(x,lambda) dgamma(y, shape, rate=1)
= k(x,shape) dgamma(y, x+shape, rate=2).
where k(x,shape) is constant as 1/{x! (shape+1)!} where rate=1.

Previously we called PrCat(x,λ) as likelihood. This time we may call dgamma(x,shape,rate) as prior distribution of λ.

Bayes’ rule

According to Bayes’ rule
Posterior = Liklihood * Prior / Evidence.
In upper case
Liklihood = LkPoi(x,λ); Pr(x|λ) Prior = dgamma(λ , shape); Pr(λ)
Evidence = ∫ PrPoi(x, λ)dλ

Posterior distribution & Maximum a posteriori estimation (MAP)

Posterior distribution of λin the world may be estimated and visualized as a graph from the data x.

#
Posterior.Poi<-function(x,a){
  Post<-NULL
  for (l in 0:max(x)){
    p<-1
    for(k in 0:max(x)){
      p<-p*dgamma(l,shape=a+k,rate=2)
    }
    Post<-rbind(Post,c(l,p))
  }
  Post<-data.frame(Post)
  names(Post)<-c("Lambda","posterior")
  sp<-sum(Post$posterior)
  Post$posterior<-Post$posterior/sp
  Post
}
#####
Post<-Posterior.Poi(x=poik,a=lambda)
plot(Post,col=3,type="l", main="Posterior distribution")
#
Post<-Posterior.Poi(x=poi20,a=lambda)
lines(Post,col=4)

If your data is not large, then the shape of the distribution will be dull also in the case of posterior estimate as in the case of the likelifood.

Posterior.Poi is not necessarily congugated. You can calculate Likelifood and prior separately and then get prior as a product of the both.

Posterior.Poi<-function(x,a){
  lam<-0:max(x)
  p<-Likelifood.Poi(x)$Likelihood*dgamma(lam,shape=a,rate=1)
  p<-p/sum(p) #1/sum(p)=Const
  Post<-data.frame(lam,p)
  names(Post)<-c("Lambda","posterior")
  Post}