如何绘制对数似然函数的图像

Question

我可以帮您翻译成中文。需要绘制-log似然函数在-pi和pi之间的图形。

-log似然函数

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)

无法绘制该函数。以下是错误消息：

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

- bowshock

---

错误信息似乎相当清晰。llh（teta，x）将仅是单个值，因为sum未矢量化。（不仅如此，您还应该查看您的高中数学书关于取负数对数的内容。） - IRTFM

上次我检查的时候，圆周率不是4。 - Spacedman

Answer 1

你的函数只适用于一个 teta 值和多个 x 值，或者多个 teta 值和一个 x 值。否则会得到错误值或警告。

例如：teta=1 和 teta=2 时的 llh：

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

不是同一个意思：

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

如果你尝试做三件事情：

> 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

这基本上源自于：

> 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

因为R试图从长度为20的向量中减去一个长度为3的向量。它会将长度为3的向量重复6次，直到完成，然后只剩下两个元素。R想：“嗯……我这样做可以，但它看起来不对，所以这里有一个警告。”

您可以重新编写该函数，使其适用于两个参数的向量，或者对该函数进行矢量化包装。

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

- Spacedman

---

感谢您的帮助。 - bowshock

Answer 2

你需要使用函数sapply（读取? sapply），因为你的代码不是向量化的。

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

将您的函数向量化，如下所示：

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

---

矢量化更有趣。 - Spacedman

向量化调用mapply，因此我会同意风格，而不是实质。 - RockScience

谢谢你的帮助。 - bowshock

Answer 3

使用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)

我同意spaced man说向量化很有趣！

- KH Kim