I want to fit a standardized Student's-t distribution. The log-likelihood is given by:

\begin{align*} log \mathcal{L}(\nu | l_1,...,l_n)=\sum_{i=1}^n \left( log \left( (\pi (\nu-2))^{-\frac{1}{2}}\Gamma \left(\frac{\nu}{2} \right)^{-1} \Gamma \left(\frac{\nu+1}{2} \right) \left(1+\frac{l^2}{\nu-2} \right)^{\text{$-\frac{1+\nu}{2}$}}\right)\right) \end{align*}

I try do to this with the optim command. My R code is (data):

# log likelihood

loglikstandardizedt <-function(par){
if(par>0) return(-sum(log((pinumber*(par-2))^(-1/2)*gamma(par/2)^(-1)*gamma((par+1)/2)*(1+standresidalvewma^2/(par-2))^(-(1+par)/2))))
else return(Inf)

optim(startvalue,loglikstandardizedt, method="BFGS")
param = optim(startvalue,loglikstandardizedt, method="BFGS")$par

I can control it via looking at the plot:

  denstiystandtresid<-function (x) (pinumber*(param-2))^(-1/2)*gamma(param/2)^(-1)*gamma((param+1)/2)*(1+x^2/(param-2))^(-(1+param)/2)


I have the problem, that the output of optim, that means the value for "par" is highly dependend on the starting value. What is the right starting value? And what is the best value for my $\nu$?

The variation is very large, I get values from 27 up to 215, so this is a huge difference?


In general, when using optim for one dimensional problem, you should set the method argument to Brent otherwise optim is indeed unreliable.

There are several problems with your code as written:

  • pinumber you can use pi instead (it's a reserved constant in R)
  • par is a reserved keyword in R so it's better to use another name.
  • loglikstandardizedt assumes a vector standresidalvewma which is not declared in the call to loglikstandardizedt
  • 1
    $\begingroup$ ok, but now it seems I have to insert a lower and upper value, how can I do this and what values should I use? $\endgroup$ – Stat Tistician Apr 24 '13 at 11:17
  • 1
    $\begingroup$ is e.g optim(90, fn=loglikstandardizedt, method="Brent",lower=4,upper=250) ok? $\endgroup$ – Stat Tistician Apr 24 '13 at 11:19
  • 1
    $\begingroup$ First of all, the time complexity of univariate optimization is sub-linear in the width of initial interval between lower and upper so there is very little cost in taking an interval that is too large. clearly the formula you use is only meaningful when ν>2 so that's your lower bound. As for the upper bound, I concur that values of ν beyond, say 100, should are for all intend and purposes yield similar results so that would be a tentative upper bound. $\endgroup$ – user603 Apr 24 '13 at 13:21
  • $\begingroup$ To the OPs defense, stats::optim even uses par as a named argument and so do other optimization functions in R like those in the optimx package. Though I agree that it would be nice to avoid. $\endgroup$ – Benjamin Christoffersen Nov 3 '18 at 11:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.