r - ggplot2 stat_function plots the wrong function -

i want plot loglikelihood function of series of independent bernoulli distributed random variables y parameter p being function (the logistic function) of feature x. logistic function has parameter b. parameter want estimate. want plot loglikelihood function of b. want in r using ggplot2, because want become better in those.

my creation of loglikelihood function can , should done better, not point here. problem plotted loglikelihood constant on interval (-5,5). seems wrong. especially, because when call function arbitrary b in interval returns different value. why happen? thank you.

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   (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 

the problem loglikelihood function. must pass "vectorized" funciton stat_function. functions in r return vector if pass in vector. example sin(1:10) return sine of numbers 1 through 10. however, when vector of values passed function, 1 value returned

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

since doesn't behave "normal" r function, have problem. easiest way fix wrap function definition in vectorize command. observe

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 ... 

now vloglikelihood behaves r function should. can plot doing

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