[R] [FORGED] How to find the likelihood, MLE and plot it?
peter dalgaard
pdalgd at gmail.com
Thu Nov 19 19:08:14 CET 2015
> 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
More information about the R-help
mailing list