[R] Fat tails

[Alberto Monteiro]

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