MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

### Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

# Thread Subject: Help on fitting a distribution to data with fminunc

 Subject: Help on fitting a distribution to data with fminunc From: Chris Date: 19 Jul, 2012 14:55:25 Message: 1 of 5 Hi All, I want to fit a skewed-t distribution to some return data (a single series of 2000 observations) in order to obtain the maximum likelihood estimates of the degrees of freedom (nu) and skewness (lambda) parameters. When I run the code below I get the following: "fminunc stopped because it cannot decrease the objective function along the current search direction. nu =   1.0e-007 *     0.1512          0 lam =     1.5154" The nu parameter is too low and I don't know if I am on track, Any suggestions on improving the code in order to get some reliable parameter estimates would be greatly appreciated. Thank you for your time! %%%%%%%%%%%%%%%%%code%%%%%%%%%%%%%%%%%%%%%%% function [nu,lambda]=EstSktParmsc(data) %% Optimisation Parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Tol = 1e-7; T = size(data,1); thita=[6;0]; % Initial estimates x0 = [thita]; %% Optimisation Options Set - fmincon options = optimset('fminunc'); options = optimset(options , 'Algorithm' , 'trust-region-reflective'); options = optimset(options , 'Display' , 'iter'); options = optimset(options , 'Diagnostics' , 'on'); options = optimset(options , 'LargeScale' , 'off'); options = optimset(options , 'MaxFunEvals' , 5000); options = optimset(options , 'MaxIter' , 5000); options = optimset(options , 'TolFun' , 1e-12); options = optimset(options , 'TolCon' , 1e-12); options = optimset(options , 'TolX' , 1e-12); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Optimisation Subroutine%%%%%%%%%%% [nu,lambda]=fminunc(@skewtdis_llf,x0,options); %% Objective function(Log was obtained from Andrew Patton's site)%%%%%%%%     function [LL, likelihoods] = skewtdis_llf(data, theta) theta=[6;0];         nu = theta(1); lambda = theta(2); [T,k] = size(data); logc = gammaln((nu+1)/2) - gammaln(nu/2) - 0.5*log(pi*(nu-2)); c = exp(logc); a = 4*lambda.*c.*((nu-2)./(nu-1)); logb = 0.5*log(1 + 3*lambda.^2 - a.^2); b = exp(logb); find1 = (data<(-a./b)); find2 = (data>=(-a./b)); LL1 = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*data(find1)+a)./(1-lambda)).^2); LL2 = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*data(find2)+a)./(1+lambda)).^2); LL = sum(LL1) + sum(LL2); LL = -LL; if nargout>1      likelihoods=zeros(size(data));      likelihoods(find1)=LL1;      likelihoods(find2)=LL2; end     end end
 Subject: Help on fitting a distribution to data with fminunc From: Matt J Date: 19 Jul, 2012 15:30:28 Message: 2 of 5 "Chris" wrote in message ... > Hi All, > > I want to fit a skewed-t distribution to some return data (a single series of 2000 observations) in order to obtain the maximum likelihood estimates of the degrees of freedom (nu) and skewness (lambda) parameters. =============== Hard to say, but here are a few issues, (1) You are calling fminunc seemingly in disregard of the input/output syntax in the documentation. In particular, lambda is one of your unknown variables, yet you are returning it as the fval output argument. All unknowns should be bundled together in a single vector. (2) The unknown loglikelihood parameters theta are passed as the 2nd argument to skewtdis_llf. FMINUNC expects it to be the first argument. A solution would be to wrap it in an anomyous function  fminunc(@(theta) skewtdis_llf(data,theta),x0,options). (3) It is not clear how the "data" variable gets passed to skewtdis_llf in the code you've shown. Similarly, your output is a variable "lam", which appears nowhere in the code. This leads me to think we're not seeing the code you're actually running, making it hard for us to help you troubleshoot... In any case, when you have problems with FMINUNC or the other solvers, you should call them with all output arguments, including exitflag, so that the output conditions can be studied.
 Subject: Help on fitting a distribution to data with fminunc From: Chris Date: 19 Jul, 2012 19:14:12 Message: 3 of 5 "Matt J" wrote in message ... > > Hard to say, but here are a few issues, > > (1) You are calling fminunc seemingly in disregard of the input/output syntax in the documentation. In particular, lambda is one of your unknown variables, yet you are returning it as the fval output argument. All unknowns should be bundled together in a single vector. > > (2) The unknown loglikelihood parameters theta are passed as the 2nd argument to skewtdis_llf. FMINUNC expects it to be the first argument. A solution would be to wrap it in an anomyous function > fminunc(@(theta) skewtdis_llf(data,theta),x0,options). > > (3) It is not clear how the "data" variable gets passed to skewtdis_llf in the code you've shown. Similarly, your output is a variable "lam", which appears nowhere in the code. This leads me to think we're not seeing the code you're actually running, making it hard for us to help you troubleshoot... > > > In any case, when you have problems with FMINUNC or the other solvers, you should call them with all output arguments, including exitflag, so that the output conditions can be studied. Dear Matt, thank you for your help, I think I have addressed the above mentioned inconsistencies and now I have to do the fine tuning. The code now has as follows: function [x,fval,exitflag]=estskewt(data) %% Objective function     function [LL, likelihoods] = skewt_llf(pars, data) nu = pars(1); %Degrees of Freedom parameter lambda = pars(2); %Skewness parameter [T,k] = size(data); logc = gammaln((nu+1)/2) - gammaln(nu/2) - 0.5*log(pi*(nu-2)); c = exp(logc); a = 4*lambda.*c.*((nu-2)./(nu-1)); logb = 0.5*log(1 + 3*lambda.^2 - a.^2); b = exp(logb); find1 = (data<(-a./b)); find2 = (data>=(-a./b)); LL1 = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*data(find1)+a)./(1-lambda)).^2); LL2 = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*data(find2)+a)./(1+lambda)).^2); LL = sum(LL1) + sum(LL2); LL = -LL; if nargout>1      likelihoods=zeros(size(data));      likelihoods(find1)=LL1;      likelihoods(find2)=LL2; end     end %% Optimisation Parameters & Starting Values nu=5; lambda=0; pars=[nu;lambda]; parsL=[0;-15]; parsU=[Inf;15]; % Initial estimates x0 = [5;0]; %% Optimisation Options Set - fmincon optopt=optimset('Algorithm', 'interior-point', 'MaxIter',100,'Display','Iter','LargeScale','on') %% Optimisation Subroutine [x,fval,exitflag]=fmincon(@skewt_llf,pars, [],[],[],[],parsL,parsU,[],optopt,data); pars=x; [LL,likelihoods]=skewt_llf(pars,data); end After running it with my data (15 iterations) what I get is: "Local minimum possible. Constraints satisfied. fmincon stopped because the size of the current step is less than the default value of the step size tolerance and constraints were satisfied to within the default value of the constraint tolerance. parameters =     2.0001    -0.0314 fval =  -5.8077e+003 exitflag =      2"
 Subject: Help on fitting a distribution to data with fminunc From: Matt J Date: 21 Jul, 2012 16:17:29 Message: 4 of 5 "Chris" wrote in message ... > > After running it with my data (15 iterations) what I get is: > "Local minimum possible. Constraints satisfied. > > fmincon stopped because the size of the current step is less than > the default value of the step size tolerance and constraints were > satisfied to within the default value of the constraint tolerance. > > > > > parameters = > > 2.0001 > -0.0314 > > > fval = > > -5.8077e+003 > > > exitflag = > > 2" ================= Ok, but what exactly is the issue now? FMINCON is telling you that it thinks it succeeded. Incidentally, by returning "likelihoods" as the 2nd output argument from your objective function, you are preventing yourself from turning on the 'GradObj' option and supplying your own calculated gradient. Indeed, it looks like computing the analytical gradient wouldn't be a hard calculation and would reduce the time FMINCON spends on finite differences.
 Subject: Help on fitting a distribution to data with fminunc From: Chris Date: 21 Jul, 2012 22:38:34 Message: 5 of 5 "Matt J" wrote in message ... > "Chris" wrote in message ... > > Incidentally, by returning "likelihoods" as the 2nd output argument from your objective function, you are preventing yourself from turning on the 'GradObj' option and supplying your own calculated gradient. Indeed, it looks like computing the analytical gradient wouldn't be a hard calculation and would reduce the time FMINCON spends on finite differences. Thanks a million, I think it is ok, parameter estimates make sense.

## Tags for this Thread

### What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us