I have a process that generates distributions that look like normal, 
but with fat tails. I would like to model them as if they were
a mix of two normals (with zero mean). Is there any function that
recovers the standard deviations and the probability based on a 
sample?
I would like to generate a sample like this:
r.mydist <- function(n, prob, mean1, sd1=1, mean2, sd2=1) {
  # aesthetical stuff
  if (missing(mean1) & missing(mean2)) mean1 <- 0
  if (missing(mean1)) mean1 <- mean2
  if (missing(mean2)) mean2 <- mean1
  # now, let's work
  # test is true if we must use the first normal (mean1, sd1)
  test <- runif(n) < prob
  # create the return value
  rval <- 0
  rval[n] <- 0
  # generate the first normal when appropriate
  if (any(test))
    rval[test] <- rnorm(sum(test), mean1, sd1)
  # generate the second normal when appropriate
  if (any(!test))
    rval[!test] <- rnorm(sum(!test), mean2, sd2)
  # game over
  return(rval)
}
# test
x <- r.mydist(1000, 0.1, 0, 1, 0, 5)
hist(x)
if (require(fBasics))
  cat("the sample kurtosis is", kurtosis(x), "\n")
Now, with this sample x, how can I get back sd1, sd2 and prob?
Alberto Monteiro