2014-06-26 86 views
4

我想繪製一系列獨立伯努利分佈隨機變量y的對數似然函數,其中參數p是一些特徵x的函數(邏輯函數)。這個邏輯函數也有一個參數b。這是我想估計的參數。所以我想繪製作爲b函數的對數似然函數。我想使用ggplot2在R中這樣做,因爲我想在這些中變得更好。ggplot2 stat_function繪製了錯誤的函數

我的loglikelihood函數的創建可以也應該做得更好,但這不是我的觀點。問題是繪製的對數似然性在區間(-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 

回答

3

問題是你的loglikelihood函數。你必須將一個「向量化」的功能傳遞給stat_function。如果你傳入一個向量,R中的大多數函數將返回一個向量。例如sin(1:10)將通過10返回數字1的正弦然而,當值的矢量傳遞給你的函數,只有一個返回值是

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") 

enter image description here