[R] [FORGED] How to find the likelihood, MLE and plot it?

[peter dalgaard]

> On 19 Nov 2015, at 18:12 , David Winsemius <dwinsemius at comcast.net> wrote:
> 
>> 
>> On Nov 19, 2015, at 7:35 AM, C W <tmrsg11 at gmail.com> wrote:
>> 
>> ah.  Let me fix that and get back to you.
>> 
>> On a side note, why can't I put fn inside curve(), why do I have to use
>> Vectorize()?
>> 
>> Here's the code and the message:
> 
>>>>>> fn <- function(theta){
>>>>>> 
>>>>>> sum(0.5 * (xvec - rep(theta, 7)) ^ 2 / 1 + 0.5 * log(1))
>>>>>> 
>>>>>> }
> 
>> 
>>> curve(fn, -5, 20)
>> Error in curve(fn, -5, 20) :
>> 'expr' did not evaluate to an object of length ’n’
> 
> 
> The sum function is not vectorized. In oder for a function to be considered “vectorized”, it needs to return a vector of the same length as its input.
> 

Even so, you would have to consider the implications of (xvec - rep(theta, 7)) if theta is a vector. You should have a fighting chance with something like

fn <- function(theta)
  colSums(outer(xvec, theta, dnorm, sd=1, log=TRUE))

or, if you are squeamish about relying on argument order:

fn <- function(theta)
colSums(outer(xvec, theta, function(x,th) dnorm(x=x, mean=th, sd=1, log=TRUE) ))

-pd

> — 
> David.
>> 
>> 
>> On Thu, Nov 19, 2015 at 10:29 AM, peter dalgaard <pdalgd at gmail.com> wrote:
>> 
>>> 
>>> On 19 Nov 2015, at 16:17 , C W <tmrsg11 at gmail.com> wrote:
>>> 
>>>> Hi Rolf,
>>>> 
>>>> I think the MLE should be 1.71, no?  And yes, I am aware of the
>>>> maximum=TRUE argument.  I still feel something is wrong here.
>>>> 
>>> 
>>> Just read more carefully what Rolf said: Your fn is MINUS the
>>> log-likelihood. So the graph is upside-down.
>>> 
>>> -pd
>>> 
>>> 
>>>> Thanks!
>>>> 
>>>> On Wed, Nov 18, 2015 at 6:23 PM, Rolf Turner <r.turner at auckland.ac.nz>
>>>> wrote:
>>>> 
>>>>> On 19/11/15 11:31, C W wrote:
>>>>> 
>>>>>> Dear R list,
>>>>>> 
>>>>>> I am trying to find the MLE of the likelihood function.  I will plot
>>> the
>>>>>> log-likelihood to check my answer.
>>>>>> 
>>>>>> Here's my R code:
>>>>>> 
>>>>>> xvec <- c(2,5,3,7,-3,-2,0)
>>>>>> 
>>>>>> fn <- function(theta){
>>>>>> 
>>>>>> sum(0.5 * (xvec - rep(theta, 7)) ^ 2 / 1 + 0.5 * log(1))
>>>>>> 
>>>>>> }
>>>>>> 
>>>>>> gn <- Vectorize(fn)
>>>>>> 
>>>>>> curve(gn, -5, 20)
>>>>>> 
>>>>>> optimize(gn, c(-5, 20))
>>>>>> 
>>>>>> $minimum
>>>>>> 
>>>>>> [1] 1.714286
>>>>>> 
>>>>>> $objective
>>>>>> 
>>>>>> [1] 39.71429
>>>>>> 
>>>>>> 
>>>>>> The MLE using optimize() is 1.71, but what curve() gives me is the
>>>>>> absolute
>>>>>> minimum.
>>>>>> 
>>>>>> I think 1.71 is the right answer, but why does the graph showing it's
>>> the
>>>>>> minimum?  What is going on here?
>>>>>> 
>>>>> 
>>>>> Your graph shows that there is indeed a *minimum* at 1.71.  And
>>> optimise()
>>>>> is correctly finding that minimum.
>>>>> 
>>>>> If you want optimise() to find the maximum, set maximum=TRUE.  In which
>>>>> case it will return "20" (or something very close to 20).
>>>>> 
>>>>> Your function fn() appears not to be the log likelihood that you had in
>>>>> mind.  Perhaps you the negative of fn()???
>>>>> 
>>>>> cheers,
>>>>> 
>>>>> Rolf Turner
>>>>> 
>>>>> --
>>>>> Technical Editor ANZJS
>>>>> Department of Statistics
>>>>> University of Auckland
>>>>> Phone: +64-9-373-7599 ext. 88276
>>>>> 
>>>> 
>>>>     [[alternative HTML version deleted]]
>>>> 
>>>> ______________________________________________
>>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>>> and provide commented, minimal, self-contained, reproducible code.
>>> 
>>> --
>>> Peter Dalgaard, Professor,
>>> Center for Statistics, Copenhagen Business School
>>> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
>>> Phone: (+45)38153501
>>> Office: A 4.23
>>> Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>> 
>> 	[[alternative HTML version deleted]]
>> 
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
> 
> David Winsemius
> Alameda, CA, USA

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com