This is something commonly done using either the Pearson or Johnson family of distributions. (It has been many years since I looked at this, so there may be other choices by now too.) The MATLAB stats TB offers the Pearson family (which I think I recall is the better choice anyway.)
lookfor pearson
pearscdf - Pearson cumulative distribution function
pearspdf - Pearson probability density function
pearsrnd - Pearson system random numbers
For example, given the first 4 moments, we can use those tools.
help pearscdf
PEARSCDF Cumulative distribution function (CDF) for Pearson system of distributions
F = PEARSCDF(X, MU, SIGMA, SKEW, KURT) returns an array F of the
Pearson system's cumulative distribution function evaluated at
double/single array X. The evaluated Pearson distribution has mean MU,
standard deviation SIGMA, skewness SKEW, and kurtosis KURT. MU and SIGMA
can be double/single scalars or arrays. If they are arrays, each must
have the same size as X. SKEW, and KURT must be double/single
scalars.
The kurtosis KURT must be greater than the square of the skewness SKEW
plus 1. The kurtosis of the normal distribution 3.
F = PEARSCDF(__, 'upper') returns the complement of the CDF. For all
Pearson types except 4, this syntax uses an algorithm that more
accurately computes the extreme upper-tail probabilities.
[F, TYPE] = PEARSCDF(__) returns the type TYPE of the specified
distribution within the Pearson system. Type is a scalar integer from
0 to 7. The eight distribution types in the Pearson system correspond
to the following distributions:
Type 0: Normal distribution
Type 1: Four-parameter beta
Type 2: Symmetric four-parameter beta
Type 3: Three-parameter gamma
Type 4: Not related to any standard distribution. Density proportional
to (1+((x-mu)/sigma)^2)^(-a) * exp(-b*arctan((x-mu)/sigma))
where a and b are quantities related to the differential
equation that defines the Pearson system.
Type 5: Inverse gamma location-scale
Type 6: F location-scale
Type 7: Student's t location-scale
[F,TYPE,COEFS] = PEARSCDF(__) returns the coefficients of the quadratic
polynomial that defines the distribution via the following differential
equation:
d(log(p(x)))/dx = -(a + x) / (COEFS(1) + COEFS(2)*x + COEFS(3)*x^2)
Example: Evaluate the Pearson distribution CDF at parameters that
correspond to the standard Normal:
% mu=0, sigma=1, skew=0, kurt=3 for the Pearson distribtion
% corresponds to the standard Normal distribution
% Evaluate this CDF from -1 to 1 in increments of 0.1
[cdf, type] = pearscdf(-1:0.1:1, 1, 2, 0, 3)
See also PEARSPDF, PEARSRND, CDF
Documentation for pearscdf
doc pearscdf
How would you use them? Simple enough. For example, I know the normal distribution has skewness zero, and kurtosis 3. So a normal distribution, with mean 1, and standard deviation 1/2 would have the cdf...
[~,T] = pearscdf([],1,1/2,0,3)
A Pearson type 0 is a normal distribution. Whoopsie do! I got it right!
And we can plot the cdf as:
fplot(@(x) pearscdf(x,1,1/2,0,3))
You can use pearsrnd to generate random numbers from that distribution, etc.