數理統計 example of EM 與 Example 7.2.13

[jackdan (Sean)]

example of EM:
Σ用R稍微算一下就可以得到答案了，不過必須在重新
scale過一次，程式附給大家參考一下
==================================================
>x=c(6, 7, 5, 6, 0, 2, 1, 1, 3, 6, 2, 5)
>T=matrix(x,byrow=T,nrow=3)
>cov(t(T))/cov(t(T))[1,3]
     [,1] [,2] [,3]
[1,] 0.50 0.25 1.00
[2,] 0.25 0.50 0.75
[3,] 1.00 0.75 2.50
==================================================




Example 7.2.13:
很遺憾的，我看完大家寫完的作業後只有兩個人真的有
用數值解的方式去找出答案，剩下的我用我的直覺猜測
你們應該都是把 Example 7.2.9 的東西重寫過一便，
然後就用一句"透過R"得知答案。我也曾經唸過碩一，
所以我能理解你們的心情。另外有一些人是寫"透過R模
擬"，這就比較誇張了，助教斗膽請教一下學弟妹，這
題目在那裡需要生成隨機變數了，並不是有要寫程式的部
份就要稱為模擬。當然有兩位同學附上了程式，不過稍
微看了一下覺得同學在程式設計的邏輯上還可以再改進
一下。

至於為何我會覺得大家都沒用寫程式去驗證呢？理由其
實是課本的答案寫錯了，然後大家就跟著一起錯了。所
以我附上我寫的程式給大家看，讓大家參考一下

=================================================
X <- c(16,18,22,25,28)

likelihood <- function(temp){
  k <- temp[1]
  p <- mean(X)/k
  1/log(prod(choose(k,X)*p^X*(1-p)^(k-X)))
}

tmp <- optim(max(X),likelihood,method="L-BFGS-B",
         lower=c(max(X),0),upper=c(Inf,1))$par

k <- floor(tmp[1])+1
p <- mean(X)/k

cat("k=",k,"\np=",p,"\n")
==================================================
最後的結果是 k=191, p=0.1141, likelihood=5.112894e-07
至於例題中的答案 k=190, p=0.1147, likelihood=5.112893e-07
算法不困難，prod(dbinom(X,k,mean(X)/k))就是likelihood
很明顯的可以看到k=191時的likelihood略大於k=190



最後附上習題 7.5 的解法：
稍微整理一下式子後可以得到方程式
-380x^3+419x^2-40x+15/16=0
有人用牛頓法去解，有人用R中的uniroot function
不過最簡單的方法是用 polyroot
程式如下
=================================================
> x=polyroot(c(15/16,-40,419,-380))
> x=as.real(x)
Warning message:
imaginary parts discarded in coercion
> k=floor(1/x)
> k
[1] 26 15  1
==================================================


最後心得：請大家記得 random variable 是大寫，observation
才是小寫。不要混淆了

--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.115.45.118
※ 編輯: jackdan         來自: 140.115.45.118       (03/08 13:18)
推 jack52014:感謝助教
推 sintheta: 感謝助教
推 kensmile26: 謝助教