ggplot2的stat_function绘制了错误的函数。

Question

我希望绘制一系列独立伯努利分布随机变量y的对数似然函数，该随机变量的参数p是某个特征x的函数（逻辑函数）。这个逻辑函数还有一个参数b。这就是我想要估计的参数。所以我想将对数似然作为b的函数来绘制。我想使用ggplot2在R中进行操作，因为我想更好地掌握它们。

我的对数似然函数的创建和完善可以做得更好，但这不是我的重点。问题在于绘制的对数似然在区间(-5,5)上是常数。这似乎是错误的。尤其是当我用该区间内的任意b值调用函数时，返回的值总是不同的。这是为什么呢？谢谢。

library(ggplot2)
set.seed(123)

# parameters 
n=100
mu=0
s=2
b<-0.2

# functions 
logit <- function(x,b){1/(1+exp(-b*x))}

# simulation of data
x<-rnorm(n,mu,s)
y_prob<-logit(x,b)
y<-rbinom(n,1,y_prob)
df<-data.frame(x,y)

# loglikelihood function 
loglikelihood<-function(b,df){
  prd<-1
  for (i in 1:NROW(df)){
    events<-logit(df$x[i],b)
    nonevents<-1-events
    prd<-prd*events^df$y[i]*nonevents^(1-df$y[i])
  }
  return(sum(log(prd)))
}


loglikelihood(0.3,df)

p2<-ggplot(data=data.frame(b=c(-5,5)), aes(b)) + stat_function(fun=loglikelihood, args=list(df=df))
p2<-p2+xlab("b") + ylab("loglikelihood")
p2

- Vincent Lous

Answer 1

问题在于你的对数似然函数。你必须向stat_function传递一个“矢量化”的函数。如果你向大多数R中的函数传递一个矢量，它们将返回一个矢量。例如，sin（1：10）将返回数字1到10的正弦值。但是，当一个值的矢量被传递到你的函数时，只有一个值会被返回。

loglikelihood(seq(-5,5, by=.1), df)
# [1] -20534.44

因为它的行为不像一个“正常”的R函数，所以你遇到了这个问题。最简单的解决方法是将你的函数定义包裹在 Vectorize 命令中。观察下面的例子：

vloglikelihood <- Vectorize(loglikelihood, vectorize.args="b")
vloglikelihood(seq(-5,5, by=.1), df)
# [1] -463.67919 -454.67142 -445.66980 -436.67470 -427.68654 -418.70574 ...

现在vloglikelihood的表现就像一个良好的R函数一样。然后我们可以像您正在做的那样将其绘制出来。

ggplot(data=data.frame(b=c(-5,5)), aes(b)) + 
    stat_function(fun=vloglikelihood, args=list(df=df)) +
    xlab("b") + ylab("loglikelihood")

- MrFlick