3

I want to graph the log likelihood function between -pi and pi.

the log likelihood function

llh <- function (teta,x) {

  sum(log((1-cos(x-teta))/(2*pi)))
}

x=c(3.91,4.85,2.28,4.06,3.70,4.04,5.46,3.53,2.28,1.96,2.53,3.88,2.22,3.47,4.82,2.46,2.99,2.54,0.52,2.50)

teta=seq(-4,4, by=0.01)

y = llh(teta,x)

plot(teta, llh(teta,x), pch=16)

It failed to plot the function. Here's the error message:

Warning message:
In x - teta :
  longer object length is not a multiple of shorter object length
>   
> plot(teta, llh(teta,x), pch=16)
Error in xy.coords(x, y, xlabel, ylabel, log) : 
  'x' and 'y' lengths differ
In addition: Warning message:
In x - teta :
  longer object length is not a multiple of shorter object length
Spacedman
  • 92,590
  • 12
  • 140
  • 224
bowshock
  • 31
  • 1
  • 1
  • 2
  • The error message seems fairly clear. `llh(teta, x)` will only be a single value because `sum` is not vectorized. (Not only that but you should review what your high school math book says about taking the log of a negative number.) – IRTFM Feb 10 '14 at 07:53
  • Last time I checked, pi wasn't 4. – Spacedman Feb 10 '14 at 08:08

3 Answers3

6

As written your function will work for one value of teta and several x values, or several values of teta and one x values. Otherwise you get an incorrect value or a warning.

Example: llh for teta=1 and teta=2:

> llh(1,x)
[1] -34.88704>
> llh(2,x)
[1] -60.00497

is not the same as:

> llh(c(1,2),x)
[1] -49.50943

And if you try and do three:

> llh(c(1,2,3),x)
[1] -49.52109
Warning message:
In x - teta :
  longer object length is not a multiple of shorter object length

which fundamentally comes from:

> cos(x-c(1,2,3))
[...]
Warning message:
In x - c(1, 2, 3) :
  longer object length is not a multiple of shorter object length

because R is trying to subtract a length-3 vector from a length-20 vector. It repeats the length-3 vector over the length-20 vector until its done it six times, then it only has two elements left. Hmmmm... thinks R, I'll do this but it looks wrong so here's a warning.

You can either rewrite the function to work on vectors for both arguments, or vectorise the function by wrapping it.

> vllh = Vectorize(llh,"teta")
> vllh(c(1,2,3),x)
[1] -34.88704 -60.00497 -67.30765
> plot(teta, vllh(teta,x))

ll plot

Spacedman
  • 92,590
  • 12
  • 140
  • 224
2

you'll need to use the function sapply (read ?sapply) as your code is not vectorised

plot(teta, sapply(X=teta, FUN=function(teta) llh(teta, x=x)), type="l")

OR

to vectorise your function as such: 

llh2 <- function (teta,x) {
  sapply(X=teta, FUN=function(teta) sum(log((1-cos(x-teta))/(2*pi))) )
}

plot(teta, llh2(teta,x), type="l")
RockScience
  • 17,932
  • 26
  • 89
  • 125
0

Another solution using outer:

llh <- function (teta,x) {
  X <- outer(x, teta, "-")
  Y <- log((1-cos(X))/(2*pi))
  apply(Y, 2, sum)
}

x=c(3.91,4.85,2.28,4.06,3.70,4.04,5.46,3.53,2.28,1.96,2.53,3.88,2.22,3.47,4.82,2.46,2.99,2.54,0.52,2.50)

teta=seq(-4,4, by=0.01)  

plot(teta, llh(teta,x), pch=16)

I agree with spaced man that Vectorize is interesting!

KH Kim
  • 1,155
  • 1
  • 7
  • 14