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)