Maximum entropy sampler

Question

I want to sample from a distribution which has fixed to a given values mean(=0), standard deviation(=1), skewness(=0) and kurtosis. I also want this distribution to be as general as possible, i.e. to have Kullback-Leibler divergence from uniform distribution as small as possible (this condition is equivalent to principle of Maximum Entropy), just like the normal distribution would be for kurtosis=3.

I know, that most probably there is no hope to have closed form for such distribution in general case. I am only interested in sampling from it. I accept reasonable numerical approximations.

I am not overly concerned with efficiency - I can wait about 2 days to have just 2000 samples.

Years age I wrote a sort of genetic algorithm to address this problem:

1. Start with random sample of 2000 values sampled from uniform distribution, which will be later called population.
2. Create many (about 100) variants of the population by first removing random subsample of - say - 200 samples from the population, and then reinserting the same number of random numbers sampled from uniform distribution.
3. Find the variant which has mean, sd, skewness and kurtosis as close as possible to the target parameters (choice of metric should not be critical to this algorithm, since all conditions are independent from each other).
4. Continue on step 2 with population being the best chosen variant from step 3.

The algorithm is slow, but in the end it yields what I understand a good approximation of maximum entropy distribution.

Is this algorithm correct? Or is there any better method for getting general leptokurtic distributions?

Answer 1

You can look for a discrete distribution, with the desired first four moments, and with the maximum entropy possible. You can then interpolate the cumulative distribution function to sample from it.

In R, it can be done as follows.

kurtosis <- 3
n <- 100
x <- seq(-5,5,length=n)
dx <- mean(diff(x))
# Opposite of the Entropy, to minimize
f <- function(p) sum( p * log(p) )
# The first moments
g <- function(p) c( sum(p)*dx, sum(x*p)*dx, sum(x^2*p)*dx, sum(x^3*p)*dx, sum(x^4*p)*dx )
# Maximize the entropy subject to those constraints
library(Rsolnp)
r <- solnp( 
  rep(1/n,n), 
  f,           # Function to minimize
  eqfun = g,   # Equality constraints
  eqB   = c(1, mean=0, var=1, skewness=0, kurtosis),
  LB=rep(0,n), UB=rep(1,n) 
)
# Beware: it is not very precise at the boundaries of the interval
plot(x, r$pars, type="l", log="y", las=1)
lines(x, dnorm(x), lty=3)
# Sample from the corresponding distribution
q <- approxfun( c(0,cumsum(r$pars)*dx), c(x[1]-dx,x) )
r <- function(n) q(runif(n))
qqnorm(r(1e4))

Answer 2

If you only have the kurtosis issue to address, you can use Student $t$-distribution with $\nu$ degrees of freedom that has kurtosis of $6/(\nu-4)$ for $\nu>4$. You would also need to normalize the variance to 1 (it is equal to $\nu/(\nu-2)$ for the original Student distribution).

If you are going to be OK with discrete distributions, you can have a distribution with support on $\{-1/2, -x, x, 1/2\}$ with corresponding probabilities $(1-p)/2, p/2, p/2, (1-p)/2$. It has 0 odd moments, its variance is $p x^2 + (1-p)/4 = 1/4 + p(x^2-1/4)$, so the relation between $x$ and $p$ would be $$p(x^2-1/4)=3/4, \quad p=\frac3{4x^2-1}, \quad x=\sqrt{\frac3{4p}+\frac14}.$$ Finally, it has the fourth moment $$1/16 (1-p)/2 + p x^4 = \frac{1-p}{32} + \frac{(3p+1)^2}{16p}=\frac{19p^2+11p+2}{32p},$$ which is also its kurtosis by virtue of the unit variance, for $0\le p\le 1$.