Bernoulli likehood function

Let us define the Bernoulli likelihood function:

bernoulliLikelihood <- function(observedSample, piGreek) {
    n <- length(observedSample)
    (1 - piGreek)^(n - sum(observedSample)) * piGreek^sum(observedSample)
}

Suppose to have the following random sample:

tenTosses <- c(1, 0, 0, 1, 0, 0, 0, 0, 0, 1)

Likelihood values for three different values of the parameter:

bernoulliLikelihood(tenTosses, 0.1)
[1] 0.0004782969
bernoulliLikelihood(tenTosses, 0.3)
[1] 0.002223566
bernoulliLikelihood(tenTosses, 0.5)
[1] 0.0009765625

It is possible to compute the three likelihood values using only a function call:

bernoulliLikelihood(tenTosses, c(0.1, 0.3, 0.5))
[1] 0.0004782969 0.0022235661 0.0009765625

To plot the likelihood function we prepare a vector of values for the parameter spanning along the parameter space:

abscissa <- seq(0, 1, 0.01)

Then we compute and plot the likelihood function:

ordinate <- bernoulliLikelihood(tenTosses, abscissa)
plot(abscissa, ordinate, type = "l")

The relative (normalized) likelihood function is easily computed dividing by the maximum:

plot(abscissa, ordinate/max(ordinate), type = "l")

The relative likelihood function is useful to compare likelihood function for different sample sizes. Let us suppose to have a sample of 100 tosses with the same relative frequency of successes:

oneHundredTosses <- c(rep(0, 70), rep(1, 30))

Using the relative likelihood function we can compare the likelihood for the two samples on the same plot:

ordinate2 <- bernoulliLikelihood(oneHundredTosses, abscissa)
plot(abscissa, ordinate/max(ordinate), type = "l")
lines(abscissa, ordinate2/max(ordinate2), col = "red")

Bernoulli log-likehood function

For computation it is convenient to use the logarith transformation of the likelihood function:

bernoulliLogLikelihood <- function(observedSample, piGreek) {
    n <- length(observedSample)
    (n - sum(observedSample)) * log(1 - piGreek) + sum(observedSample) * log(piGreek)
}

Here is the log-likelihood plot:

logOrdinate <- bernoulliLogLikelihood(tenTosses, abscissa)
plot(abscissa, logOrdinate, type = "l")

And here is the relative log-likelihood:

plot(abscissa, logOrdinate - max(logOrdinate), type = "l")

As for the likelihood function, the relative log-likelihood allows to easily compare the shape of the function for different sample sizes:

logOrdinate2 <- bernoulliLogLikelihood(oneHundredTosses, abscissa)
plot(abscissa, logOrdinate - max(logOrdinate), type = "l")
lines(abscissa, logOrdinate2 - max(logOrdinate2), col = "red")

LS0tCnRpdGxlOiAiTGlrZWhvb2QgYW5kIGxvZy1saWtlbGlob29kIGZ1bmN0aW9ucyBmb3IgdGhlIEJlcm5vdWxsaSBtb2RlbCIKYXV0aG9yOiAiRG9tZW5pY28gVmlzdG9jY28iCmRhdGU6ICJCUyBDbGFzcyAtIExlY3R1cmUgb2YgMTctMTgvMTAvMjAxNiIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyMgQmVybm91bGxpIGxpa2Vob29kIGZ1bmN0aW9uCkxldCB1cyBkZWZpbmUgdGhlIEJlcm5vdWxsaSBsaWtlbGlob29kIGZ1bmN0aW9uOgoKYGBge3J9CmJlcm5vdWxsaUxpa2VsaWhvb2QgPC0gZnVuY3Rpb24ob2JzZXJ2ZWRTYW1wbGUsIHBpR3JlZWspIHsKICAgIG4gPC0gbGVuZ3RoKG9ic2VydmVkU2FtcGxlKQogICAgKDEgLSBwaUdyZWVrKV4obiAtIHN1bShvYnNlcnZlZFNhbXBsZSkpICogcGlHcmVla15zdW0ob2JzZXJ2ZWRTYW1wbGUpCn0KYGBgClN1cHBvc2UgdG8gaGF2ZSB0aGUgZm9sbG93aW5nIHJhbmRvbSBzYW1wbGU6IApgYGB7cn0KdGVuVG9zc2VzIDwtIGMoMSwgMCwgMCwgMSwgMCwgMCwgMCwgMCwgMCwgMSkKYGBgCkxpa2VsaWhvb2QgdmFsdWVzIGZvciB0aHJlZSBkaWZmZXJlbnQgdmFsdWVzIG9mIHRoZSBwYXJhbWV0ZXI6CmBgYHtyfQpiZXJub3VsbGlMaWtlbGlob29kKHRlblRvc3NlcywgMC4xKQpiZXJub3VsbGlMaWtlbGlob29kKHRlblRvc3NlcywgMC4zKQpiZXJub3VsbGlMaWtlbGlob29kKHRlblRvc3NlcywgMC41KQpgYGAKSXQgaXMgcG9zc2libGUgdG8gY29tcHV0ZSB0aGUgdGhyZWUgbGlrZWxpaG9vZCB2YWx1ZXMgdXNpbmcgb25seSBhIGZ1bmN0aW9uIGNhbGw6CmBgYHtyfQpiZXJub3VsbGlMaWtlbGlob29kKHRlblRvc3NlcywgYygwLjEsIDAuMywgMC41KSkKYGBgCgpUbyBwbG90IHRoZSBsaWtlbGlob29kIGZ1bmN0aW9uIHdlIHByZXBhcmUgYSB2ZWN0b3Igb2YgdmFsdWVzIGZvciB0aGUgcGFyYW1ldGVyIHNwYW5uaW5nIGFsb25nIHRoZSBwYXJhbWV0ZXIgc3BhY2U6CmBgYHtyfQphYnNjaXNzYSA8LSBzZXEoMCwgMSwgMC4wMSkKYGBgCgpUaGVuIHdlIGNvbXB1dGUgYW5kIHBsb3QgdGhlIGxpa2VsaWhvb2QgZnVuY3Rpb246CmBgYHtyfQpvcmRpbmF0ZSA8LSBiZXJub3VsbGlMaWtlbGlob29kKHRlblRvc3NlcywgYWJzY2lzc2EpCnBsb3QoYWJzY2lzc2EsIG9yZGluYXRlLCB0eXBlID0gImwiKQpgYGAKCgpUaGUgcmVsYXRpdmUgKG5vcm1hbGl6ZWQpIGxpa2VsaWhvb2QgZnVuY3Rpb24gaXMgZWFzaWx5IGNvbXB1dGVkIGRpdmlkaW5nIGJ5IHRoZSBtYXhpbXVtOgpgYGB7cn0KcGxvdChhYnNjaXNzYSwgb3JkaW5hdGUvbWF4KG9yZGluYXRlKSwgdHlwZSA9ICJsIikKYGBgCgpUaGUgcmVsYXRpdmUgbGlrZWxpaG9vZCBmdW5jdGlvbiBpcyB1c2VmdWwgdG8gY29tcGFyZSBsaWtlbGlob29kIGZ1bmN0aW9uIGZvciBkaWZmZXJlbnQgc2FtcGxlIHNpemVzLiBMZXQgdXMgc3VwcG9zZSB0byBoYXZlIGEgc2FtcGxlIG9mIDEwMCB0b3NzZXMgd2l0aCB0aGUgc2FtZSByZWxhdGl2ZSBmcmVxdWVuY3kgb2Ygc3VjY2Vzc2VzOgpgYGB7cn0Kb25lSHVuZHJlZFRvc3NlcyA8LSBjKHJlcCgwLCA3MCksIHJlcCgxLCAzMCkpCmBgYAoKVXNpbmcgdGhlIHJlbGF0aXZlIGxpa2VsaWhvb2QgZnVuY3Rpb24gd2UgY2FuIGNvbXBhcmUgdGhlIGxpa2VsaWhvb2QgZm9yIHRoZSB0d28gc2FtcGxlcyBvbiB0aGUgc2FtZSBwbG90OgpgYGB7cn0Kb3JkaW5hdGUyIDwtIGJlcm5vdWxsaUxpa2VsaWhvb2Qob25lSHVuZHJlZFRvc3NlcywgYWJzY2lzc2EpCnBsb3QoYWJzY2lzc2EsIG9yZGluYXRlL21heChvcmRpbmF0ZSksIHR5cGUgPSAibCIpCmxpbmVzKGFic2Npc3NhLCBvcmRpbmF0ZTIvbWF4KG9yZGluYXRlMiksIGNvbCA9ICJyZWQiKQpgYGAKCiMjIEJlcm5vdWxsaSBsb2ctbGlrZWhvb2QgZnVuY3Rpb24KRm9yIGNvbXB1dGF0aW9uIGl0IGlzIGNvbnZlbmllbnQgdG8gdXNlIHRoZSBsb2dhcml0aCB0cmFuc2Zvcm1hdGlvbiBvZiB0aGUgbGlrZWxpaG9vZCBmdW5jdGlvbjoKYGBge3J9CmJlcm5vdWxsaUxvZ0xpa2VsaWhvb2QgPC0gZnVuY3Rpb24ob2JzZXJ2ZWRTYW1wbGUsIHBpR3JlZWspIHsKICAgIG4gPC0gbGVuZ3RoKG9ic2VydmVkU2FtcGxlKQogICAgKG4gLSBzdW0ob2JzZXJ2ZWRTYW1wbGUpKSAqIGxvZygxIC0gcGlHcmVlaykgKyBzdW0ob2JzZXJ2ZWRTYW1wbGUpICogbG9nKHBpR3JlZWspCn0KYGBgCgpIZXJlIGlzIHRoZSBsb2ctbGlrZWxpaG9vZCBwbG90OgpgYGB7cn0KbG9nT3JkaW5hdGUgPC0gYmVybm91bGxpTG9nTGlrZWxpaG9vZCh0ZW5Ub3NzZXMsIGFic2Npc3NhKQpwbG90KGFic2Npc3NhLCBsb2dPcmRpbmF0ZSwgdHlwZSA9ICJsIikKYGBgCgpBbmQgaGVyZSBpcyB0aGUgcmVsYXRpdmUgbG9nLWxpa2VsaWhvb2Q6CmBgYHtyfQpwbG90KGFic2Npc3NhLCBsb2dPcmRpbmF0ZSAtIG1heChsb2dPcmRpbmF0ZSksIHR5cGUgPSAibCIpCmBgYAoKQXMgZm9yIHRoZSBsaWtlbGlob29kIGZ1bmN0aW9uLCB0aGUgcmVsYXRpdmUgbG9nLWxpa2VsaWhvb2QgYWxsb3dzIHRvIGVhc2lseSBjb21wYXJlIHRoZSBzaGFwZSBvZiB0aGUgZnVuY3Rpb24gZm9yIGRpZmZlcmVudCBzYW1wbGUgc2l6ZXM6CmBgYHtyfQpsb2dPcmRpbmF0ZTIgPC0gYmVybm91bGxpTG9nTGlrZWxpaG9vZChvbmVIdW5kcmVkVG9zc2VzLCBhYnNjaXNzYSkKcGxvdChhYnNjaXNzYSwgbG9nT3JkaW5hdGUgLSBtYXgobG9nT3JkaW5hdGUpLCB0eXBlID0gImwiKQpsaW5lcyhhYnNjaXNzYSwgbG9nT3JkaW5hdGUyIC0gbWF4KGxvZ09yZGluYXRlMiksIGNvbCA9ICJyZWQiKQpgYGAKCg==