sandk<-function(x){
##
#Calculates the mean, variance, skew, and kurtosis#
# for a data set and returns them in that order.#
#The formulas for skew and kurtosis are from page 85#
# of Kendall and Steward 1969, vol.1#
##
n <- length(x)
m1p <- mean(x)
m2 <- sum((x-m1p)^2)/n
m3 <- sum((x-m1p)^3)/n
m4 <- sum((x-m1p)^4)/n
k1 <- m1p
k2 <- n*m2/(n-1)
k3 <- ((n^2)/((n-1)*(n-2)))*m3
k4 <- (((n^2)/((n-1)*(n-2)))/(n-3))*((n+1)*m4 - 3*(n-1)*m2^2)
g1 <- k3/(k2^(3/2))
g2 <- k4/(k2^2)
return(c(k1,k2,g1,g2))
}


fleishtarget<-function(x,a){
##
#The target function for solving equations 18, 19, and 20#
# from page 523 of Fleishman.#
#It does this by changing the system of three equations into a #
# minimization problem.  Set the equations equal to zero, square,#
# and sum up.  The b, c, and d that minimize the set of equations # 
# at zero must also solve the three individually.#
##
b<-x[1]
cc<-x[2]
d<-x[3]
g1<-a[1]
g2<-a[2]
(2 - ( 2*b^2 + 12*b*d + g1^2/(b^2+24*b*d+105*d^2+2)^2 + 30*d^2 ) )^2 +
(g2 - ( 24*(b*d+cc^2*(1+b^2+28*b*d)+d^2*(12+48*b*d+141*cc^2+225*d^2)) ) )^2+
(cc - (g1/(2*(b^2+24*b*d+105*d^2+2)) ) )^2
}


findbcd<-function(skew,kurtosis){
##
#Uses the built in minimization function to solve for b, c, and d#
# if the skew and kurtosis are given. Try findbcd(1.75,3.75) and#
# compare to Table 1 on page 524 of Fleishman.
##
optim(c(1,0,0),fleishtarget,a=c(skew,kurtosis),method="BFGS",
control=list(ndeps=rep(1e-10,3),reltol=1e-10,maxit=1e8))
}


rfleish<-function(n,mean=0,var=1,skew=0,kurtosis=0){
##
#Generates n random variables with specified first four moments#
# using Fleishman's power method. Note that not all combinations#
# of skew and kurtosis are possible (see Figure 1 on page 527).#  
#Must satisfy skew^2 < 0.0629576*kurtosis + 0.0717247.#
##
Z<-rnorm(n,0,1)
bcd<-findbcd(skew,kurtosis)$par
b<-bcd[1]
cc<-bcd[2]
d<-bcd[3]
a<--1*cc
Y<-a+b*Z+cc*Z^2+d*Z^3
X<-mean+sqrt(var)*Y
return(X)}