STA 410/2102, Fall 2015, Script for Assignment #1.

> source("ass1-funs.r")
> 
> options(digits=17)

Set the sample sizes and observed counts for the data in the assignment handout.

> n <- 130
> m1 <- 25
> m2 <- 25
> 
> x  <- 75
> x1 <- 20
> x2 <-  6

Create a contour plot of the log likelihood function, for reference when discussing the results.

> log_likelihood_contour <- function (n,m1,m2,x,x1,x2, ...)
+ {
+     grid <- seq(0.01, 0.99, by=0.01)
+     ll <- matrix(nrow=length(grid),ncol=length(grid))
+     for (i in 1:length(grid))
+         for (j in 1:length(grid))
+             ll[i,j] <- log_likelihood(c(grid[i],grid[j]),n,m1,m2,x,x1,x2)
+     contour (grid, grid, ll, ...)
+ }
> 
> log_likelihood_contour(n,m1,m2,x,x1,x2,nlevels=200)

plot of chunk unnamed-chunk-4

Try the methods with various initial values, to see how sensitive they are to initial values. Does only a few iterations for Newton and method of scoring, just to see if they are converging or diverging.

> init_list <- list (c(x1/m1,x2/m2),
+                    c(0.5,0.5),
+                    c(0.1,0.1),
+                    c(0.05,0.95))
> 
> for (init in init_list) {
+     cat ("\nTRYING INITIAL VALUE",init,"\n\n")
+     cat("\nAlternating maximization\n\n")
+     p_alt <- mle_alt (n,m1,m2,x,x1,x2,init)
+     cat("\nMultivariate Newton iteration\n")
+     p_mvn <- mle_mvn (n,m1,m2,x,x1,x2,init,7)
+     cat("\nMethod of scoring\n")
+     p_mos <- mle_mos (n,m1,m2,x,x1,x2,init,7)
+     cat("\nUsing nlm\n\n")
+     p_nlm <- mle_nlm (n,m1,m2,x,x1,x2,init)
+     print(p_nlm)
+ }


TRYING INITIAL VALUE 0.80000000000000004 0.23999999999999999 


Alternating maximization

i = 1 p1 = 0.84615973276776923   p2 = 0.27469089013085279 
i = 2 p1 = 0.83341349646038121   p2 = 0.28149745269102089 
i = 3 p1 = 0.83077805762026391   p2 = 0.28291450844583144 
i = 4 p1 = 0.83022377614011722   p2 = 0.28321295545116498 
i = 5 p1 = 0.83010679171269075   p2 = 0.28327596277724615 
i = 6 p1 = 0.83008208329244559   p2 = 0.28328927144453253 
i = 7 p1 = 0.83007686378814238   p2 = 0.28329208285633756 
i = 8 p1 = 0.83007576116331827   p2 = 0.28329267677114889 
i = 9 p1 = 0.8300755282312533   p2 = 0.28329280223711006 
i = 10 p1 = 0.83007547902373968   p2 = 0.28329282874212869 
i = 11 p1 = 0.83007546862851922   p2 = 0.28329283434138597 
i = 12 p1 = 0.83007546643250052   p2 = 0.28329283552424422 
i = 13 p1 = 0.83007546596858561   p2 = 0.28329283577412634 
i = 14 p1 = 0.83007546587058245   p2 = 0.28329283582691445 
i = 15 p1 = 0.83007546584987901   p2 = 0.28329283583806608 
i = 16 p1 = 0.8300754658455054   p2 = 0.28329283584042197 
i = 17 p1 = 0.83007546584458147   p2 = 0.28329283584091958 
i = 18 p1 = 0.8300754658443863   p2 = 0.28329283584102466 
i = 19 p1 = 0.830075465844345   p2 = 0.28329283584104697 
i = 20 p1 = 0.83007546584433611   p2 = 0.28329283584105169 
i = 21 p1 = 0.83007546584433434   p2 = 0.28329283584105264 
i = 22 p1 = 0.83007546584433411   p2 = 0.2832928358410528 
i = 23 p1 = 0.83007546584433389   p2 = 0.28329283584105291 
i = 24 p1 = 0.83007546584433389   p2 = 0.28329283584105291 

Multivariate Newton iteration

At iteration 1 :

x = 0.80000000000000004 0.23999999999999999 
f(x) = 14.823717948717935 14.823717948717942 

At iteration 2 :

x = 0.83428607068977545 0.27908612058634408 
f(x) = -0.86976160891783039 0.47342823773695031 

At iteration 3 :

x = 0.83015086405183069 0.28323651182752291 
f(x) = -0.017738990011661571 0.003806334172718806 

At iteration 4 :

x = 0.83007548857122737 0.28329282085591129 
f(x) = -5.6052169519205108e-06 6.6454701297402607e-07 

At iteration 5 :

x = 0.83007546584433589 0.28329283584105158 
f(x) = -4.9382720135326963e-13 6.0396132539608516e-14 

At iteration 6 :

x = 0.83007546584433389 0.28329283584105291 
f(x) = 3.5527136788005009e-15 3.5527136788005009e-15 

At iteration 7 :

x = 0.83007546584433389 0.28329283584105291 
f(x) = 3.5527136788005009e-15 3.5527136788005009e-15 

Method of scoring

At iteration 1 :

x = 0.80000000000000004 0.23999999999999999 
f(x) = 14.823717948717935 14.823717948717942 

At iteration 2 :

x = 0.83408567480423768 0.27885766927683092 
f(x) = -0.77161236262138999 0.55533650956522607 

At iteration 3 :

x = 0.82976324803289614 0.28325032333335248 
f(x) = 0.10930969254761891 0.051029116387525164 

At iteration 4 :

x = 0.83010794430802637 0.28327245194968109 
f(x) = -0.0081458386823385354 0.00070004080751928655 

At iteration 5 :

x = 0.83007262893246225 0.28329344974897508 
f(x) = 0.00086362013892227196 0.00022142662407276248 

At iteration 6 :

x = 0.83007573922780453 0.2832927112886568 
f(x) = -7.4693360510025286e-05 -5.4989056863519181e-06 

At iteration 7 :

x = 0.83007544093492192 0.28329284386753745 
f(x) = 7.2391330618870597e-06 1.3057199481636417e-06 

Using nlm

Warning: NaNs produced```

Warning: NaNs produced“

Warning: NaNs produced```