Driving test: finding the likelihood โปรแกรม R จะตอบคำถามนี้ ได้หรือไม่

Source:

ครั้งที่ 1 สอบผ่าน 150 คน (สอบผ่านในครั้งแรก)
ครั้งที่ 2 สอบผ่าน 125 คน (เคยสอบ 1 ครั้ง และไม่ผ่าน)
ครั้งที่ 3 สอบผ่าน 100 คน (เคยสอบ 2 ครั้ง และไม่ผ่าน)
ครั้งที่ 4 สอบผ่าน 75 คน (เคยสอบ 3 ครั้ง และไม่ผ่าน)
ครั้งที่ 5 และมากกว่า สอบผ่าน 50 คน (เคยสอบ 4 ครั้ง และไม่ผ่าน)

ต้องการหาค่า Maximum Likelihood Estimate จากสอบผ่าน 450 ครั้ง ไม่ผ่าน 750 ครั้ง
โปรแกรม R จะตอบคำถามนี้ ได้หรือไม่ ?

Likelihood Estimate เป็นการคาดคะเน Proportion
ของ Population จาก Sample ทีมีอยู่
สอบผ่าน ไม่ผ่าน เป็น Binomial Distribution
ถ้า N มีค่ามาก ใช้วิธี Likelihood Estimate จะมี Error
ถ้า Likelihood มีค่าเข้าใกล้ 0 เช่น 0.00000005
ตัวเลขทศนิยมส่วนท้ายๆ จะถูกตัดออก
และ Computer จะเข้าใจว่า Likelihood = 0

แก้ไขโดย take log
ถ้า likelihood = 1, log likelihood = 0
ถ้า likelihood น้อยกว่า 1, log likelihood ค่าติดลบ
เช่น ln(0.000001) = - 13.81551

(1) function log.likelihood(sequence, p)
(2) sequence, p.parameter, possible.p
(3) Plot Graph
Log Likelihood as a Function of P

กำหนดให้ P = 0.000 เพิ่มค่าครั้งละ 0.001 ถึง P = 1.000
possible.p จำนวน = 1,001 ค่า
คำนวณค่า Log Likelihood จาก sequence และทุกๆค่าของ possible.p
เช่น sequence = 1 จำนวน 450 ครั้ง, sequence = 0 จำนวน 750 คร้ง
rbinom(n, size, prob) แต่ละครั้ง ได้จำนวน 0 และ 1 ต่างกันไป

ถ้า P = 0, Log Likelihood = - 5,500
ถ้า P = 1, Log Likelihood = - 5,500
ที่ Log Likelihood = - 600 (จุดสูงสุด) P = 0.38

copy คำสั่งในตอนล่าง วางใน R console ของ โปรแกรม R
ซึ่งจะสร้าง function log.likelihood(sequence, p)
กำหนดค่า sequence, possible.p
และ plot graph "log likelihood as a function of P"
เริ่มตั้งแต่ ## ถึง ## วางใน notepad
และ copy ระหว่าง ## ถึง ## วางใน โปรแกรม R

##
log.likelihood <- function(sequence, p)
{
log.likelihood <- 0
for (i in 1:length(sequence))
{
if (sequence[i] == 1)
{
log.likelihood <- log.likelihood + log(p)
}
else
{
log.likelihood <- log.likelihood + log(1 - p)
}
}
return(log.likelihood)
}

p.parameter <- 450/1200
sequence <- rbinom(1200, 1, p.parameter)
sequence

possible.p <- seq(0, 1, by = 0.001)
jpeg('Likelihood_Concavity.jpg')
possible.p

library('ggplot2')
win.graph(width = 10, height = 10, pointsize = 10)
qplot(possible.p,
sapply(possible.p, function (p) {log.likelihood(sequence, p)}),
geom = 'line',
main = 'Log Likelihood as a Function of P',
xlab = 'P',
ylab = 'Log Likelihood')
##

Source: https://www.r-bloggers.com/doing-maximum-likelihood-estimation-by-hand-in-r/


อธิบายคำสั่งที่ใช้ใน โปรแกรม R
https://dl.dropboxusercontent.com/u/1999671/EPI_56/r/mle2/women.htm

dropbox stop render html page, Oct 2016!!
อธิบายคำสั่งที่ใช้ใน โปรแกรม R ย้ายไปที่
http://epistat.usite.pro/r/mle2/women.htm