Hello stranger, and welcome! 👋😊
I'm Rasmus Bååth, data scientist, engineering manager, father, husband, tinkerer,
tweaker, coffee brewer, tea steeper, and, occasionally, publisher of stuff I find
interesting down below👇
Update:
pingr is now on CRAN! Due to requests from the CRAN mantainers the ping() function had to be renamed beep() in order to not be confused with the Unix tool
ping. Thus the package is now called beepr instead and can be downloaded from CRAN by running the following in an R session:
install.packages("beepr")
So just replace pingr and ping with beepr and beep below, otherwise everything is the same :)
Original post:
pingr is an R package that contains one function, ping(), with one purpose: To go ping on whatever platform you are on (thanks to the
audio package). It is intended to be useful, for example, if you are running a long analysis in the background and want to know when it is ready. It’s also useful if you want to irritate colleagues. You could, for example, use ping() to get notified when your package updates have finished:
update.packages(ask=FALSE); ping()
The binomial test is arguably the conceptually simplest of all statistical tests: It has only one parameter and an easy to understand distribution for the data. When introducing null hypothesis significance testing it is puzzling that the binomial test is not the first example of a test but sometimes is introduced long after the t-test and the ANOVA (as
here) and sometimes is not introduced at all (as
here and
here). When introducing a new concept, why not start with simplest example?
It is not like there is a problem with students understanding the concept of null hypothesis significance testing too well. I’m not doing the same misstake so here follows the Bayesian First Aid alternative to the binomial test!
So I have a secret project. Come closer. I’m developing an R package that implements Bayesian alternatives to the most commonly used statistical tests. Yes you heard me, soon your t.testing days might be over! The package aims at being as easy as possible to pick up and use, especially if you are already used to the classical .test functions. The main gimmick is that the Bayesian alternatives will have the same calling semantics as the corresponding classical test functions save for the addition of bayes. to the beginning of the function name. Going from a classical test to the Bayesian version will be as easy as going from t.test(x1, x2, paired=T) to bayes.t.test(x1, x2, paired=T).
The package does not aim at being
some general framework for Bayesian inference or
a comprehensive collection of Bayesian models. This package should be seen more as a quick fix; a first aid for people who want to try out the Bayesian alternative. That is why I call the package Bayesian First Aid.
I’m playing blog ping-pong with John Kruschke’s
Doing Bayesian Data Analysis blog as he was partly inspired by my silly
post on Bayesian mascots when writing a nice piece on
Icons for the essence of Bayesian and frequentist data analysis. That piece, in turn, inspired me resulting in the following wobbly backhand.
The confidence interval is, for me, one of the more tricky statistical concepts. As opposed to p-values and posterior distributions which can be explained to a newbie “pretty easily” (especially using
John’s nice plots), I find that the concept of confidence intervals is really hard to communicate. Perhaps the easiest way to explain confidence intervals is as a bootstrap procedure (which is beautifully visualized
here). However, when I’ve tried this explanation it seems like people often end up believing that confidence intervals are more or less like
credible intervals.
Another way of explaining confidence intervals is as the region of possible null hypotheses resulting in corresponding significance tests that are not rejected. This explanation is harder to make than the bootstrap explanation, perhaps due to the double negation in the last sentence. Turns out it wasn’t easy to make a corresponding nice explanatory animation either, but that’s what I tried to do anyway…
The following animation shows the construction of a confidence interval for an ordinary least squares regression. First the the data is plotted and the least squares line is calculated. Then a null hypothesis line is randomly drawn and using this line a test is made whether this null hypothesis can be rejected or not. The scatter of thinner red lines shows lines drawn from the sampling distribution of the null hypothesis line. This is repeated a number of times and the null hypothesis lines that are not rejected form the confidence band of the regression line.
Why would Bayesian statistics need a mascot/symbol/logo? Well, why not? I don’t know of any other branch of statistics that has a mascot but many programming languages have.
R has an “R”,
Python has a python snake, and
Go has an adorable little gopher. While Bayesian statistics isn’t a programming language it could be characterized as a technique/philosophy/framework with an active community and as such it wouldn’t be completely off if it had a mascot.
However, it does not matter if you think it would be silly if Bayesian statistics had a mascot. There is already a mascot and he is not a particularly good one. Meet Mr Bayes:
So R is awesome. Felt good to get that out of the way! But sometimes I long for some small syntax additions while programming away in R in the night…
The picture is of a restaurant called Risotto in Berlin.
You’ve probably heard about the t distribution. One good use for this distribution is as an alternative to the normal distribution that is more robust against outliers. But where does the t distribution come from? One intuitive characterization of the t is as
a mixture of normal distributions. More specifically, as a mixture of an infinite number of normal distributions with a common mean $\mu$ but with precisions (the inverse of the variance) that are randomly distributed according to
a gamma distribution. If you have a hard time picturing an infinite number of normal distributions you could also think of a t distribution as a normal distribution with a standard deviation that “jumps around”.
Using this characterization of the t distribution we could generate random samples $y$ from a t distribution with a mean $\mu$, a scale $s$ and a degrees of freedom $\nu$ as:
$$y \sim \text{Normal}(\mu, \sigma) $$
$$ 1/\sigma^2 \sim \text{Gamma}(\text{shape}= \nu / 2, \text{rate} = s^2 \cdot \nu / 2)$$
This brings me to the actual purpose of this post, to show off a nifty visualization of how the t can be seen as a mixture of normals. The animation below was created by drawing 6000 samples of $1/\sigma^2$ from a $\text{Gamma}(\text{shape}= 2 / 2, \text{rate} = 3^2 \cdot 2 / 2)$ distribution and using these to construct 6000 normal distribution with $\mu=0$. Drawing a sample from each of these distributions should then be the same as sampling from a $\text{t}(\mu=0,s=3,\nu=2)$ distribution. But is it? Look for yourself:
In the
last post I showed how to use Laplace approximation to quickly (but dirtily) approximate the posterior distribution of a Bayesian model coded in R. This is just a short follow up where I show how to use
importance sampling as an easy method to shape up the Laplace approximation in order to approximate the true posterior much better. But first, what is importance sampling?
Importance Sampling
Importance sampling is a method that can be used to estimate the integral of a function even though we only have access to the density of that function at any given point. Using a further resampling step, this method can be used to produce samples that are (approximately) distributed according to the density of the given function. In Bayesian statistics this is what we often want to do with a function that describes the density of the posterior distribution. In the case when there is more than one free parameter the posterior will be multidimensional and there is nothing stopping importance sampling from working in many dimensions. Below I will only visualize the one dimensional case, however.
The first step in the importance sampling scheme is to define a proposal distribution that we believe covers most of our target distribution, the distribution that we actually want to estimate. Ideally we would want to choose a proposal distribution that is the same as the target distribution, but if actually could do this we would already be able to sample from the target distribution and importance sampling would be unnecessary. Instead we can make an educated guess. A sloppy guess is to use a uniform proposal distribution that we believe covers most of the density of the target distribution. For example, assuming the target distribution has most of its density in $[-4,4]$ we could setup the importance sampling as:
Here, the target distribution is actually a $\text{Normal}(\mu=1,\sigma=1)$ distribution, but most often we do not know its shape.
Thank you for tuning in! In this post, a continuation of
Three Ways to Run Bayesian Models in R, I will:
- Handwave an explanation of the Laplace Approximation, a fast and (hopefully not too) dirty method to approximate the posterior of a Bayesian model.
- Show that it is super easy to do Laplace approximation in R, basically four lines of code.
- Put together a small function that makes it even easier, if you just want this, scroll down to the bottom of the post.
But first a picture of
the man himself:
I argued
recently that a good way of communicating statistical models graphically was by using the convention devised by John K. Kruschke in his book
Doing Bayesian Data Analysis. John Kruschke describes these diagrams in more detail on his blog:
here,
here and
here. While I believe these kinds of diagrams are great in many ways there is a problem in that they are quite tricky to make. That is, until now!
I have put together a template for the free and open source
LibreOffice Draw which makes it simple to construct Kruschke style diagrams such as the one below:
