Maximum Likelihood Estimates and the EM Algorithms I

1
Maximum Likelihood Estimates and the EM
Algorithms I

• Henry Horng-Shing Lu
• Institute of Statistics
• National Chiao Tung University
• hslu_at_stat.nctu.edu.tw

2
Part 1Computation Tools
3
Computation Tools

• R (http//www.r-project.org/) good for
  statistical computing
• C/C good for fast computation and large data
  sets
• More http//www.stat.nctu.edu.tw/subhtml/source/t
  eachers/hslu/course/statcomp/links.htm

4
The R Project

• R is a free software environment for statistical
  computing and graphics. It compiles and runs on a
  wide variety of UNIX platforms, Windows and
  MacOS.
• Similar to the commercial software of Splus.
• C/C, Fortran and other codes can be linked and
  called at run time.
• More http//www.r-project.org/

5
Download R from http//www.r-project.org/
6
Choose one Mirror Site of R
7
Choose the OS System
8
Select the Base of R
9
Download the Setup Program
10
Install R
Double click R-icon to install R
11
Execute R
Interactive command window
12
Download Add-on Packages
13
Choose a Mirror Site
Choose a mirror site close to you
14
Select One Package to Download
Choose one package to download, like rgl or
adimpro.
15
Load Packages

• There are two methods to load packages

Method 1 Click from the menu bar
Method 2 Type library(rgl) in the command
window
16
Help in R (1)

• What is the loaded library?
• help(rgl)

17
Help in R (2)

• How to search functions for key words?
• help.search(key words)
• It will show all functions has the key words.
• help.search(3D plot)

18
Help in R (3)

• How to find the illustration of function?
• ?function name
• It will show the usage, arguments, author,
  reference, related functions, and examples.
• ?plot3d

19
R Operators (1)

• Mathematic operators
• , -, , /,
• Mod
• sqrt, exp, log, log10, sin, cos, tan,

20
R Operators (2)

• Other operators
• sequence operator
• matrix algebra
• lt, gt, lt, gt inequality
• , ! comparison
• , , , and, or
• formulas
• lt-, assignment

21
Algebra, Operators and Functions

• gt 12
• 1 3
• gt 1gt2
• 1 FALSE
• gt 1gt2 2gt1
• 1 TRUE
• gt A 13
• gt A
• 1 1 2 3
• gt A6
• 1 6 12 18
• gt A/10
• 1 0.1 0.2 0.3
• gt A2
• 1 1 0 1

gt B 46 gt AB 1 4 10 18 gt t(A)B 1 1
32 gt At(B) 1 2 3 1 4 5 6 2
8 10 12 3 12 15 18 gt sqrt(A) 1 1.000
1.1414 1.7320 gt log(A) 1 0.000 0.6931 1.0986
gt round(sqrt(A), 2) 1 1.00 1.14
1.73 gt ceiling(sqrt(A)) 1 1 2
2 gt floor(sqrt(A)) 1 1 1 1 gt eigen(At(B)) va
lues 1 3.20e01 8.44e-16 -4.09e-16 vectors
1 2 3 1, -0.2673 0.3112
-0.2353 2, -0.5345 -0.8218 -0.6637 3, -0.8018
0.4773 0.7100
22
Variable Types
23
Define Your Own Function (1)

• Use "fix(myfunction)"
• a window will show up
• function(parameter)
• statements
• return (object)
• if you want to return some values
• Save the document
• Use "myfunction(parameter)" in R

24
Define Your Own Function (2)

• Example Find all the factors of an integer

25
Define Your Own Function (3)
When you leave the program, remember to save the
work space for the next use, or the function you
defined will disappear after you close R project.
26
Read and Write Files

• Write Data to a TXT File
• Write Data to a CSV File
• Read TXT and CSV Files
• Demo

27
Write Data to a TXT File

• Usage
• write(x, file, )
• gt X matrix(16, 2, 3)
• gt X
• ,1 ,2 ,3
• 1, 1 3 5
• 2, 2 4 6
• gt write(t(X), file "d/out1.txt", ncolumns
  3)
• gt write(X, file "d/out2.txt", ncolumns 3)

d/out1.txt 1 3 5 2 4 6
d/out2.txt 1 2 3 4 5 6
28
Write Data to a CSV File

• Usage
• write.table(x, file "foo.csv", )
• gt X matrix(16, 2, 3)
• gt X
• ,1 ,2 ,3
• 1, 1 3 5
• 2, 2 4 6
• gt write.table(t(X), file "d/out1.csv", sep
  ",", col.names FALSE, row.names FALSE)
• gt write.table(X, file "d/out2.csv", sep
  ",", col.names FALSE, row.names FALSE)

d/out1.csv 1,2 3,4 5,6
d/out2.csv 1,3,5 2,4,6
29
Read TXT and CSV Files

• Usage
• read.table(file, ...)
• gt X read.table(file "d/out1.txt")
• gt X
• V1 V2 V3
• 1 1 3 5
• 2 2 4 6
• gt Y read.table(file "d/out1.csv", sep
  ",", header FALSE)
• gt Y
• V1 V2
• 1 1 2
• 2 3 4
• 3 5 6

30
Demo (1)

• Practice for read file and basic analysis
• gt Data read.table(file "d/01.csv", header
  TRUE, sep ",")
• gt Data
• Y X1 X2
• 1, 2.651680 13.808990 26.75896
• 2, 1.875039 17.734520 37.89857
• 3, 1.523964 19.891030 26.03624
• 4, 2.984314 15.574260 30.21754
• 5, 10.423090 9.293612 28.91459
• 6, 0.840065 8.830160 30.38578
• 7, 8.126936 9.615875 32.69579

01.csv
31
Demo (2)

• Practice for read file and basic analysis
• gt mean(DataY)
• 1 4.060727
• gt boxplot(DataY)
• gt boxplot(Data)

32
Part 2Motivation Examples
33
Example 1 in Genetics (1)

• Two linked loci with alleles A and a, and B and b
• A, B dominant
• a, b recessive
• A double heterozygote AaBb will produce gametes
  of four types AB, Ab, aB, ab

34
Example 1 in Genetics (2)

• Probabilities for genotypes in gametes

35
Example 1 in Genetics (3)

• Fisher, R. A. and Balmukand, B. (1928). The
  estimation of linkage from the offspring of
  selfed heterozygotes. Journal of Genetics, 20,
  7992.
• More
• http//en.wikipedia.org/wiki/Genetics
  http//www2.isye.gatech.edu/brani/isyebayes/bank/
  handout12.pdf

36
Example 1 in Genetics (4)
37
Example 1 in Genetics (5)

• Four distinct phenotypes
• AB, Ab, aB and ab.
• A the dominant phenotype from (Aa, AA, aA).
• a the recessive phenotype from aa.
• B the dominant phenotype from (Bb, BB, bB).
• b the recessive phenotype from bb.
• AB 9 gametic combinations.
• Ab 3 gametic combinations.
• aB 3 gametic combinations.
• ab 1 gametic combination.
• Total 16 combinations.

38
Example 1 in Genetics (6)

• Let , then

39
Example 1 in Genetics (7)

• Hence, the random sample of n from the offspring
  of selfed heterozygotes will follow a multinomial
  distribution
• We know that
• and
• So

40
Example 1 in Genetics (8)

• Suppose that we observe the data of
• which is a random sample from
• Then the probability mass function is

41
Estimation Methods

• Frequentist Approaches
• http//en.wikipedia.org/wiki/Frequency_probabilit
  y
• Method of Moments Estimate (MME)
• http//en.wikipedia.org/wiki/Method_of_moments_2
  8statistics29
• Maximum Likelihood Estimate (MLE)
• http//en.wikipedia.org/wiki/Maximum_likelihood
• Bayesian Approaches
• http//en.wikipedia.org/wiki/Bayesian_probability
  

42
Method of Moments Estimate (MME)

• Solve the equations when population moments are
  equal to sample moments
• for k 1, 2, , t, where t is the
  number of parameters to be estimated.
• MME is simple.
• Under regular conditions, the MME is consistent!
• More http//en.wikipedia.org/wiki/Method_of_momen
  ts_28statistics29

43
MME for Example 1

• Note MME cant assure

44
MME by R

• gt MME lt- function(y1, y2, y3, y4)
• n y1y2y3y4
• phi1 4.0(y1/n-0.5)
• phi2 1-4y2/n
• phi3 1-4y3/n
• phi4 4.0y4/n
• phi (phi1phi2phi3phi4)/4.0
• print("By MME method")
• return(phi) print(phi)
• gt MME(125, 18, 20, 24)
• 1 "By MME method"
• 1 0.5935829

45
MME by C/C
46
Maximum Likelihood Estimate (MLE)

• Likelihood
• Maximize likelihood Solve the score equations,
  which are setting the first derivates of
  likelihood to be zeros.
• Under regular conditions, the MLE is consistent,
  asymptotic efficient and normal!
• More
• http//en.wikipedia.org/wiki/Maximum_likelihood

47
Example 2 (1)

• We toss an unfair coin 3 times and the random
  variable is
• If p is the probability of tossing head, then

48
Example 2 (2)

• The distribution of of tossing head

49
Example 2 (3)

• Suppose we observe the toss of 1 heads and 2
  tails, the likelihood function becomes
• One way to maximize this likelihood function is
  by solving the score equation, which sets the
  first derivative to be zero

50
Example 2 (4)

• The solution of p for the score equation is 1/3
  or 1.
• One can check that p1/3 is the maximum point.
  (How?)
• Hence, the MLE of p is 1/3 for this example.

51
MLE for Example 1 (1)

• Likelihood
• MLE

52
MLE for Example 1 (2)
A
B
C
53
MLE for Example 1 (3)

• Checking
• 1.
• 2.
• 3. Compare ?

54
Use R to find MLE (1)

• gt MLE
• gt y1 125 y2 18 y3 20 y4 24
• gt f lt- function(phi)
• ((2.0phi)/4.0)y1 ((1.0-phi)/4.0)(y2y3)
  (phi/4.0)y4
• gt plot(f, 1/4, 1, xlab expression(varphi), ylab
  "likelihood function multipling a constant")
• gt optimize(f, interval c(1/4, 1), maximum T)
• maximum
• 1 0.5778734
• objective
• 1 7.46944e-82

55
Use R to find MLE (2)
56
Use C/C to find MLE (1)
57
Use C/C to find MLE (2)
58
Exercises

• Write your own programs for those examples
  presented in this talk.
• Write programs for those examples mentioned at
  the following web page
• http//en.wikipedia.org/wiki/Maximum_likelihood
  
• Write programs for the other examples that you
  know.

59
More Exercises (1)

• Example 3 in genetics
• The observed data are
• where , , and fall in
• such that
• Find the likelihood function and score equations
  for , , and .

60
More Exercises (2)

• Example 4 in the positron emission tomography
  (PET)
• The observed data are
• and
• The values of are known and the unknown
  parameters are .
• Find the likelihood function and score equations
  for .

61
More Exercises (3)

• Example 5 in the normal mixture
• The observed data are random
  samples from the following probability density
  function
• Find the likelihood function and score equations
  for the following parameters