# Maximum entropy sampler

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?

-
You need to define the class of distributions of interest -- maximum entropy among the distributions that are absolutely continuous on the real line? with support on a finite range? with support on the finite number of points (as your solution goes)? –  StasK Jan 17 '13 at 14:36

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

-
Thank you for such elaborate answer. The code works fine, but at the first try the quality of the fit looks poor. I'll evaluate it more and post the results when I get them. –  Adam Ryczkowski Apr 18 '12 at 11:55
You need to increase $n$ drastically, say to the order of $10^4 - 10^6$, so as to get a reasonably accurate distribution. I was going to recommend the same thing under a different name (empirical likelihood), so +1 is warranted here –  StasK May 17 '12 at 16:57
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$.
I never said it had a maximum entropy, and I really doubt that it does. I said that if the only problem you want to address is kurtosis, this could be a way to go; lower moments are easy to fix. Given how poorly your problem is defined, and how crude your initial solution is (simulation with $n=100$), I don't see why you are concerned with a very specific method to solve it (maximum entropy). –  StasK Jan 17 '13 at 14:35