I have a function set up to calculate the likelihood of a distribution. I am getting negative values and I am wondering if my function is right? and if so how can I plot my function as a curve?
y <- function(mu, x) {
n <- length(x)
-n*mu + log(mu)*sum(x) - sum(lfactorial(x))
}
x<-c(5,4,0,1,4,7,3,5,4)
mu<-c(3,4,2.5)
y(x,mu)
Any advice would be very helpful.
Below you can find the full expression of the log-likelihood from a Poisson distribution. Additionally, I simulated data from a Poisson distribution using rpois to test with a mu equal to 5, and then recover it from the data optimizing the loglikelihood using optimize
#set seed
set.seed(777)
#loglikeliood of poisson
log_like_poissson <- function(y) {
n <- length(y)
function(mu) {
log(mu) * sum(y) - n * mu - sum(lfactorial(y))
}
}
# Data simulation: Poisson with lambda = 5
y <- rpois(n=10000, lambda = 5)
# Optimization of the loglikelihood
optimise(log_like_poissson(y),
interval = c(0, 100),
maximum = TRUE)
#$maximum
#[1] 4.994493
#
#$objective
#[1] -22033.2
This code is highly based on Chapter 10 of Advanced R where you can find an extensive discussion about how to optimize the likelihood described above.
[EDITED]
For the graph part of your question, you can use the following code to see how your loglikelihood behaves at different values of mu. As you can see from the graph, the maximum of the function is at the value of mu equal of 5 (as expected).
library(ggplot2)
values_for_mu<- seq(from=0.05, to = 10 , by =0.05 )
#new loglikelihood (only depends on mu)
log_like_poissson2 <- function(mu) {
n <- length(y)
(log(mu) * sum(y)) - (n * mu) - sum(lfactorial(y))
}
#Evaluate the loglikelihood at different values of mu
values_log_like <- unlist(lapply(values_for_mu,
FUN = log_like_poissson2))
#generate a dataframe to ggplot2
df <- data.frame(values_for_mu, values_log_like)
# Plot
ggplot(df, aes(x=values_for_mu, y=values_log_like)) +
geom_line() + geom_vline(xintercept = 5, linetype="dotted",
color = "red", size=1.5) +
xlab("mu") + ylab("Value of Log-likelihood")
I don't know if your function is right, because I don't know what you want. However, I did see that your result is a sum of two negative numbers and one positive (unless that lfactorial() does something special; I don't know what it is). However, I think I can help you with the curve.
I don't know if you're familiar with the package ggplot2, but I learned a lot about it here. Anyway, how to do the line curve would be described in a section of that website, at least if you want to literally connect the dots. If you want a smooth curve, that it's this place.
Anyway, your code would be:
library(ggplot2)
yourData <- y(x,mu)
ggplot()+
geom_line(aes(x = x, y = yourData))
That's the basics, which would output:
By the way, I'm guessing that x vector is your x values, which would go in the x axis, specially because the other way was a mess and this was kind of pretty. If it's not, you can change that.
With a little more customising, you could do:
ggplot()+
geom_line(aes(x = x, y = yourData), colour = "blue", size = 1.5)+
theme_minimal()+
xlab("Whatever this is")+
ylab("What are your results?")+
ggtitle("This is your graph. Enjoy!")
PS: I don't quite understand your function, but you seem to, so maybe these graphs help you visualize your outputs and see if they look how they're supposed to.
