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)