# Create distribution model from samples?

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?

• Are you aware of R? link – tharen Jan 26 '12 at 0:43
• aware, yes. Haven't used it though. Any suggestions along that line? – DeepSpace101 Jan 26 '12 at 23:44

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)

• Actually I wanted to create the distribution model, THEN perform a monte carlo simulation on the same model to extend it out. I guess the followup of going from a kernel + random number between [-1,1] -> random sample from distribution is a separate question altogether! – DeepSpace101 Jan 26 '12 at 23:50
• The logspline package has tools to randomly generate data from a distribution fit using the 'logspline' function. Generating a random value from a kernal density estimate is as simple as randomly choosing a value from your original data, then generating a random value from the kernal used centered around the chosen datapoint. The logspline approach is a bit simpler, but neither are complicated. – Greg Snow Jan 27 '12 at 15:41
• Thanks Greg, upvoted! This is great because it solves my underlying problem in a different manner. Since I'm still new to this so you have an R code-snippet that demonstrates your solution? – DeepSpace101 Jan 27 '12 at 18:12
• I added some example code for using logspline in the answer above. – Greg Snow Jan 27 '12 at 20:19

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.

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