@@ -36,9 +36,9 @@ ind <- 0
3636
3737for (i in 1 : length(n )) {
3838 logL <- 0
39- for (j in 1 : m ) {
40- logL <- logL + lgamma(n [i ]+ 1 )- lgamma(k [j ]+ 1 )- lgamma(n [i ]- k [j ]+ 1 )
41- logL <- logL + k [j ]* log(theta [i ])+ (n [i ]- k [j ])* log(1 - theta [i ])
39+ for (j in 1 : data $ m ) {
40+ logL <- logL + lgamma(n [i ]+ 1 )- lgamma(data $ k [j ]+ 1 )- lgamma(n [i ]- data $ k [j ]+ 1 )
41+ logL <- logL + data $ k [j ]* log(theta [i ])+ (n [i ]- data $ k [j ])* log(1 - theta [i ])
4242 }
4343 if (logL > cc ) {
4444 ind <- i
@@ -52,24 +52,23 @@ layout(matrix(c(2,0,1,3),2,2,byrow=T), width=c(2/3, 1/3), heights=c(1/3, 2/3))
5252xhist <- hist(n , plot = F )
5353yhist <- hist(theta , plot = F )
5454top <- max(c(xhist $ counts , yhist $ counts ))
55- xrange <- c(0 , nmax )
55+ xrange <- c(0 , data $ nmax )
5656yrange <- c(0 , 1 )
5757
58- par(mar = c(5 ,5 , 1 , 1 ))
59- plot(n ,theta ,xlim = xrange , ylim = yrange ,ylab = " " , xlab = " " )
58+ par(mar = c(5 , 5 , 1 , 1 ))
59+ plot(n , theta , xlim = xrange , ylim = yrange ,ylab = " " , xlab = " " )
6060axis(1 )
6161mtext(" Number of Surveys" , side = 1 ,line = 2.25 , cex = 1.2 )
6262axis(2 , cex = 1.2 )
6363las = 0
6464mtext(" Rate of Return" , side = 2 ,line = 2.25 , cex = 1.2 )
6565las = 1
66- points(mean(n ),mean(theta ), col = " red" , lwd = 3 , pch = 4 ) # expectation
67- points(n [ind ],theta [ind ], col = " green" , lwd = 3 , pch = 10 ) # Maximum Likelihood
66+ points(mean(n ), mean(theta ), col = " red" , lwd = 3 , pch = 4 ) # expectation
67+ points(n [ind ], theta [ind ], col = " green" , lwd = 3 , pch = 10 ) # Maximum Likelihood
6868
69- par(mar = c(0 ,4 , 1 , 1 ))
69+ par(mar = c(0 , 4 , 1 , 1 ))
7070barplot(xhist $ counts , axes = FALSE , ylim = c(0 , top ), space = 0 ,col = " lightblue" )
7171
72- par(mar = c(4 ,0 , 1 , 3 ))
72+ par(mar = c(4 , 0 , 1 , 3 ))
7373barplot(yhist $ counts , axes = FALSE , xlim = c(0 , top ), space = 0 , horiz = TRUE ,
7474 col = " lightblue" )
75-
0 commit comments