Fitting t-distribution in R: scaling parameter

How do I fit the parameters of a t-distribution, i.e. the parameters corresponding to the 'mean' and 'standard deviation' of a normal distribution. I assume they are called 'mean' and 'scaling/degrees of freedom' for a t-distribution?

The following code often results in 'optimization failed' errors.

library(MASS)
fitdistr(x, "t")


Do I have to scale x first or convert into probabilities? How best to do that?

• It fails not because you have to scale parameters, but because optimizer fails. See my answer below. Feb 10 '16 at 15:26

fitdistr uses maximum-likelihood and optimization techniques to find parameters of a given distribution. Sometimes, especially for t-distribution, as @user12719 noticed, the optimization in the form:

fitdistr(x, "t")


fails with an error.

In this case you should give optimizer a hand by providing starting point and lower bound to start searching for optimal parameters:

fitdistr(x, "t", start = list(m=mean(x),s=sd(x), df=3), lower=c(-1, 0.001,1))


Note, df=3 is your best guess at what an "optimal" df could be. After providing this additional info your error will be gone.

Couple of excerpts to help you better understand the inner mechanics of fitdistr:

For the Normal, log-Normal, geometric, exponential and Poisson distributions the closed-form MLEs (and exact standard errors) are used, and start should not be supplied.

...

For the following named distributions, reasonable starting values will be computed if start is omitted or only partially specified: "cauchy", "gamma", "logistic", "negative binomial" (parametrized by mu and size), "t" and "weibull". Note that these starting values may not be good enough if the fit is poor: in particular they are not resistant to outliers unless the fitted distribution is long-tailed.

• Both answers (Flom and Bushmanov) are helpful. I'm picking this one, because it makes it more explicit that with the right initial values and constraints 'fitdistr' optimization converges. Feb 10 '16 at 23:01

MASS, the book (4th edition, page 110) advises against trying to estimate $$\nu$$, the degrees of freedom parameter in the $$t$$-distribution with maximum likelihood (with some literature references: Lange et al. (1989), "Robust statistical modeling Using the t distribution", JASA, 84, 408, and Fernandez & Steel (1999), "Multivariate Student-t regression models: Pitfalls and inference", Biometrika, 86, 1).

The reason is that the likelihood function for $$\nu$$ based on the t density function, may be unbounded and will in those cases not give a well defined maximum. Let us look at an artificial example where location and scale is known (as the standard $$t$$-distribution) and only the degrees of freedom is unknown. Below is some R code, simulating some data, defining the log-likelihood function and plotting it:

set.seed(1234)
n <- 10
x <- rt(n,  df=2.5)

make_loglik  <-  function(x)
Vectorize( function(nu) sum(dt(x, df=nu,  log=TRUE)) )

loglik  <-  make_loglik(x)
plot(loglik,  from=1,  to=100,  main="loglikelihood function for df     parameter", xlab="degrees of freedom")
abline(v=2.5,  col="red2")


If you play around with this code, you can find some cases where there is a well-defined maximum, especially when the sample size $$n$$ is large. But is the maximum likelihood estimator then any good?

Let us try some simulations:

t_nu_mle  <-  function(x) {
loglik  <-  make_loglik(x)
res  <-  optimize(loglik, interval=c(0.01, 200), maximum=TRUE)\$maximum
res
}

nus  <-  replicate(1000, {x <- rt(10, df=2.5)
t_nu_mle(x) }, simplify=TRUE)

> mean(nus)
[1] 45.20767
> sd(nus)
[1] 78.77813


Showing the estimation is very unstable (looking at the histogram, a sizable portion of the estimated values is at the upper limit given to optimize of 200).

Repeating with a larger sample size:

nus  <-  replicate(1000, {x <- rt(50, df=2.5)
t_nu_mle(x) }, simplify=TRUE)
> mean(nus)
[1] 4.342724
> sd(nus)
[1] 14.40137


which is much better, but the mean is still way above the true value of 2.5.

Then remember that this is a simplified version of the real problem where location and scale parameters also have to be estimated.

If the reason of using the $$t$$-distribution is to "robustify", then estimating $$\nu$$ from the data may well destroy the robustness.

• Your conclusion that the problems of estimating df might actually work against the reason to choose a t-distribution in the first place (i.e. robustness) is thought provoking. Feb 12 '16 at 20:00
• (+1) "Unbounded above" isn't a wrong answer & might well be useful for some purposes when coupled with an interval estimate. The important thing is not to blindly use the observed Fisher information to form Wald confidence intervals. May 22 '16 at 11:40

In the help for fitdistr is this example:

fitdistr(x2, "t", df = 9)


indicating that you just need a value for df. But that assumes standardization.

For more control, they also show

mydt <- function(x, m, s, df) dt((x-m)/s, df)/s
fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, lower = c(-Inf, 0))


where the parameters would be m = mean, s = standard deviation, df = degrees of freedom

• I guess I am confused about the parameters of a t-distribution. Does it have 2 (mean,df) or 3 (mean, standard deviation, df) parameters? I was wondering if one could fit the parameter 'df'. Dec 12 '12 at 21:10
• @user12719 The Student's-t distribution has three parameters: location, scale and degrees of freedom. They are not referred as mean, standard deviation and df because the mean and the variance of this distribution depend on the three parameters. Also, they do not exists in some cases. Peter Flom is fixing the df but this can be considered as an unknown parameter as well.
– user10525
Dec 12 '12 at 21:25
• @PeterFlom In the case of the Cauchy distribution it is explicit that m and s are the location and scale. I agree the notation m and s suggests that they represent the mean and standard deviation, respectively. But this may just be a simplification of \mu and \sigma as well. +1 long ago, by the way.
– user10525
Dec 12 '12 at 22:09
• @PeterFlom Does this citation from the R's help file imply that df is always 9 for students distribution? Don't you think df should be estimated as well? Actually, the absence of df is the cause for the error, and the right answer should provide some recipe to finding it. Feb 9 '16 at 17:23
• @PeterFlom BTW, if you read the help file couple of lines above your citation, you'll find why df=9 is good in their example and irrelevant here. Feb 9 '16 at 17:30

You can use fitdistrplus library after extending the location and scaling parameters for the student t in base R according to this article on wikipedia. Below is sample code

library(fitdistrplus)
x<-rt(100,23)
dt_ls <- function(x, df=1, mu=0, sigma=1) 1/sigma * dt((x - mu)/sigma, df)
pt_ls <- function(q, df=1, mu=0, sigma=1)  pt((q - mu)/sigma, df)
qt_ls <- function(p, df=1, mu=0, sigma=1)  qt(p, df)*sigma + mu
rt_ls <- function(n, df=1, mu=0, sigma=1)  rt(n,df)*sigma + mu
fit.t<-fitdist(x, 't_ls', start =list(df=1,mu=mean(x),sigma=sd(x)))
summary(fit.t)