Create distribution model from samples?

Question

If I have a set of samples, say 100-200 samples and I'd like to create the distribution model from this list of samples, what is a reasonably efficient way of doing it? Are there any opensource / easily accessible statistical libraries?

eg: if I assume it's a normal distribution, I can easily find the mean and variance and call it a day. However, the distribution might actually have non-negligible 3rd (skewness) or even 4th moment (Kurtosis) so my "assumption" is not very accurate. My gut feeling is this may work:

// assume samples already in samples[]
float avg = CalculateAverage(samples);
float variance = CalculateVariance(samples);
float skewness = CalculateSkewness(samples);
float kurtosis = CalculateKurtosis(samples);
string definedBy = "average, variance";

if(skewness > skewThreshold)
   definedBy += ", skewness";
if(kurtosis > kurtosisThreshold)
   definedBy += ", kurtosis";
Console.Writeline("This distribution is defined by " + definedBy);

So:

Q1: Will it achieve the purpose of classifying the distribution/creating the distribution model ?

Q2: What values of skewThreshold and kurtosisThreshold are reasonable?

Answer 1

If you want to decide on a specific family of distributions then you should really consider the science that leads to your data.

If you just want an estimation or feel for what the distribution looks like then you can use tools like kernal density estimation or logspline density estimation (or others). Both of the above can be done using R which is free and open source. The density function will do the kernal density estimation and the logspline function in the package of the same name will do the log spline estimation.

Here is some sample R code for the logspline method (this assumes that you have already installed the logspline package):

library(logspline)

par(mfrow=c(3,1))
hist(iris$Petal.Width)

myfit <- logspline( iris$Petal.Width, lbound=0 )
plot(myfit)

new.Petal.Width <- rlogspline(1000, myfit)
hist(new.Petal.Width)

Answer 2

You can model your distribution as a mixture of gaussians ("normal mixture"): for instance, the mixture of two gaussian distributions with the same mean and different variances has a higher kurtosis; likewise, the mixture of two gaussian distributions with different means has a non-zero skewness. If you model your distribution as a mixture of two gaussians, you need four parameters, which could be your first four moments.

In R, you can use the flexmix package to model mixtures of gaussians. If you want statistical tests to check if the skewness and kurtosis are significantly different from those of a gaussian distribution (the "thresholds" you mention), you can check the moments package.

Answer 3

I recently discovered the fitdistrplus library, so I thought I'd work up an example. The example uses some random data scaled to [-1,1] as you suggested and then naively assumes a beta distribution to generate 1000 random samples. I added a plot to illustrate the results. I use a beta distribution here, but R has many others that may work for you. The fitdist object returns parameter estimates for the requested distribution function.

I was going to suggest bootstrap resampling. However, I'm not sure how it would fit with your end objectives.

#using the fitdistrplus package from CRAN
library(fitdistrplus)

#using fake data from a beta distribution here...
#generate some skewed random numbers scaled to [-1,1]
sample.data = rbeta(n=30,shape1=5,shape2=3) * 2 - 1

#estimate the best fit assuming a beta distribution [0,1]
sample.scaled = (sample.data + 1.0) / 2.0
fit = fitdist(sample.scaled,'beta', method='mle')
summary(print(fit))

##Generate a diagnostic plot of the fit
#plot(fit)

#generate a density curve for the normal distribution
fit.x = seq(0.0,1.0,length=40)
a = fit$estimate[1]
b = fit$estimate[2]
fit.pdf = dbeta(fit.x,shape1=a,shape2=b)/2

#generate a random sample from the fitted distribution scaling to [-1,1]
fit.random = rbeta(1000,shape1=a,shape2=b)*2-1

#display the original data, PDF, and sample data
h=hist(sample.data, xlim=c(-1,1), breaks=20, freq=F
       , main='Samples from Distribution'
       ,xlab = '')
lines(fit.x*2-1,fit.pdf, col='blue',lty=2)
lines(density(fit.random),col='red')

#add a legend
leg.txt = c('sample.data','fit.pdf','fit.random')
leg.cols = c('black','blue','red')
leg.lty = c(1,2,1)
legend(-1,max(h$intensities),leg.txt, col=leg.cols, lty=leg.lty)