{"id":487,"date":"2013-07-31T21:35:56","date_gmt":"2013-08-01T01:35:56","guid":{"rendered":"http:\/\/www.clayford.net\/statistics\/?p=487"},"modified":"2014-05-17T09:36:54","modified_gmt":"2014-05-17T13:36:54","slug":"a-probability-problem-in-heredity-part-3","status":"publish","type":"post","link":"https:\/\/www.clayford.net\/statistics\/a-probability-problem-in-heredity-part-3\/","title":{"rendered":"A Probability Problem in Heredity \u2013 Part 3"},"content":{"rendered":"

In my previous two posts I showed worked solutions to problems 2.5<\/a> and 11.7<\/a> in Bulmer’s Principles of Statistics<\/a><\/em>, both of which involve the characteristics of self-fertilizing hybrid sweet peas. It turns out that problem 11.8 also<\/em> involves this same topic, so why not work it as well for completeness. The problem asks us to assume that we were unable to find an explicit solution for the maximum likelihood equation in problem 11.7 and to solve it by using the following iterative method:<\/p>\n

\\( \\theta_{1} = \\theta_{0} + \\frac{S(\\theta_{0})}{I(\\theta_{0})} \\)<\/p>\n

where \\( S(\\theta_{0}) \\) is the value of \\( \\frac{d \\log L}{d\\theta}\\) evaluated at \\( \\theta_{0}\\) and \\( I(\\theta_{0})\\) is the value of \\( -E(\\frac{d^{2}\\log L}{d\\theta^{2}})\\) evaluated at \\( \\theta_{0}\\).<\/p>\n

So we begin with \\( \\theta_{0}\\) and the iterative method returns \\( \\theta_{1}\\). Now we run the iterative method again starting with \\( \\theta_{1}\\) and get \\( \\theta_{2}\\):<\/p>\n

\\( \\theta_{2} = \\theta_{1} + \\frac{S(\\theta_{1})}{I(\\theta_{1})} \\)<\/p>\n

We repeat this process until we converge upon a value. This is called the Newton-Raphson method<\/a>. Naturally this is something we would like to have the computer do for us.<\/p>\n

First, recall our formulas from problem 11.7:<\/p>\n

\\( \\frac{d \\log L}{d\\theta} = \\frac{1528}{2 + \\theta} – \\frac{223}{1 – \\theta} + \\frac{381}{\\theta} \\)
\n\\( \\frac{d^{2}\\log L}{d \\theta^{2}} = -\\frac{1528}{(2 + \\theta)^{2}} -\\frac{223}{(1 – \\theta)^{2}} -\\frac{381}{\\theta^{2}} \\)<\/p>\n

Let’s write functions for those in R:<\/p>\n

\r\n# maximum likelihood score\r\nmls <- function(x) {\r\n\t1528\/(2 + x) - 223\/(1 - x) + 381\/x\r\n\t}\r\n# the information\r\ninf <- function(x) {\r\n\t-1528\/((2 + x)^2) - 223\/((1 - x)^2) - 381\/(x^2)\r\n\t}\r\n\r\n<\/pre>\n
Now we can use those functions in another function that will run the iterative method starting at a trial value:<\/p>\n
\r\n# newton-raphson using expected information matrix\r\nnr <- function(th) {\r\n prev <- th\r\n repeat {\r\n   new <- prev + mls(prev)\/-inf(prev)\r\n   if(abs(prev - new)\/abs(new) <0.0001)\r\n     break\r\n   prev <- new\r\n  }\r\nnew\r\n}\t\r\n<\/pre>\n
This function first takes its argument and names it \"prev\". Then it starts a repeating loop. The first thing the loop does it calculate the new value using the iterative formula. It then checks to see if the difference between the new and previous value - divided by the new value - is less than 0.0001. If it is, the loop breaks and the \"new\" value is printed to the console. If not, the loop repeats. Notice that each iteration is hopefully converging on a value. As it converges, the difference between the \"prev\" and \"new\" value will get smaller and smaller. So small that dividing the difference by the \"new\" value (or \"prev\" value for that matter) will begin to approach 0.<\/p>\n
To run this function, we simply call it from the console. Let's start with a value of \\( \\theta_{0} = \\frac{1}{4}\\), as the problem suggests:<\/p>\n
\r\nnr(1\/4)\r\n[1] 0.7844304\r\n<\/pre>\n
There you go! We could make the function tell us a little more by outputting the iterative values and number of iterations. Here's a super quick and dirty way to do that:<\/p>\n
\r\n# newton-raphson using expected information matrix\r\nnr <- function(th) {\r\n k <- 1 # number of iterations\r\n v <- c() # iterative values\r\n  prev <- th\r\n  repeat {\r\n    new <- prev + mls(prev)\/-inf(prev)\r\n    v[k] <- new\r\n    if(abs(prev - new)\/abs(new) <0.0001)\r\n     break\r\n    prev <- new\r\n    k <- k + 1\r\n    }\r\nprint(new) # the value we converged on\r\nprint(v) # the iterative values\r\nprint(k) # number of iterations\r\n}\r\n<\/pre>\n
Now when we run the function we get this:<\/p>\n
\r\nnr(1\/4)\r\n[1] 0.7844304\r\n[1] 0.5304977 0.8557780 0.8062570 0.7863259 0.7844441 0.7844304\r\n[1] 6\r\n<\/pre>\n
We see it took 6 iterations to converge. And with that I think I've had my fill of heredity problems for a while.<\/p>\n","protected":false},"excerpt":{"rendered":"
In my previous two posts I showed worked solutions to problems 2.5 and 11.7 in Bulmer’s Principles of Statistics, both… Read more<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[3,7,13],"tags":[25,28],"class_list":["post-487","post","type-post","status-publish","format-standard","hentry","category-mle","category-probability","category-using-r","tag-bulmer","tag-newton-raphson"],"_links":{"self":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/487","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/comments?post=487"}],"version-history":[{"count":10,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/487\/revisions"}],"predecessor-version":[{"id":649,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/487\/revisions\/649"}],"wp:attachment":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/media?parent=487"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/categories?post=487"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/tags?post=487"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}