Bioops

Bioinformatics=(ACGAAG->AK)+(#!/bin/sh)+(P(A|B)=P(B|A)*P(A)/P(B))

Estimate Gamma Distribution Parmaters Using MME and MLE

| Comments

This post shows how to estimate gamma distribution parameters using (a) moment of estimation (MME) and (b) maximum likelihood estimate (MLE).

The probability density function of Gamma distribution is

The MME:

We can calculate the MLE of $ \alpha $ using the Newton-Raphson method.

For $ k =1,2,…,$

where

Use the MME for the initial value of $ \alpha^{(0)} $, and stop the approximation when $ \vert \hat{\alpha}^{(k)}-\hat{\alpha}^{(k-1)} \vert < 0.0000001 $. The MLE of $ \beta $ can be found by $ \hat{\beta} = \bar{X} / \hat{\alpha} $.

Below is the R code.

gamma.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

# (a) MME
gamma_MME<-function(x){
  n<-length(x)
  mean_x<-mean(x)
  alpha<-n*(mean_x^2)/sum((x-mean_x)^2)
  beta<-sum((x-mean_x)^2)/n/mean_x
  estimate_MME<-data.frame(alpha,beta)
  return(estimate_MME)
}


# (b) MLE
gamma_MLE<-function(x){
  n<-length(x)
  mean_x<-mean(x)

  # initiate the convergence and alpha value
  converg<-1000
  alpha_prev<-gamma_MME(x)$alpha

  # initiate two vectors to store alpha and beta in each step
  alpha_est<-alpha_prev
  beta_est<-mean_x/alpha_prev

  # Newton-Raphson
  while(converg>0.0000001){
    #first derivative of alpha_k-1
    der1<-n*log(alpha_prev/mean_x)-n*digamma(alpha_prev)+sum(log(x))
    #second derivative of alpha_k-1
    der2<-n/alpha_prev-n*trigamma(alpha_prev)
    #calculate next alpha
    alpha_next<-alpha_prev-der1/der2
    # get the convergence value
    converg<-abs(alpha_next-alpha_prev)
    # store estimators in each step
    alpha_est<-c(alpha_est, alpha_next)
    beta_est<-c(beta_est, mean_x/alpha_next)
    # go to next alpha
    alpha_prev<-alpha_next
  }

  alpha<-alpha_next
  beta<-mean_x/alpha_next
  estimate_MLE<-data.frame(alpha,beta)

  return(estimate_MLE)
}

# apply
x<-rgamma(100,2,scale=5)
gammma_MME(x)
gamma_MLE(x)

Comments