Parameters seem to not be optimizing in fmincon

I have code that is using fmincon to find parameters of interest, but one parameter updates over time (in Bayesian-esque fashion) as well as in a general parameter vector. I currently am having the data generated to test the function, and will eventually input observed data (the issues below happen regardless). The code runs and a minimum is found, but shockinfo and tau end up really large and equal to each other. Additonally as Tau updates, its value equals some variant of 1/n. Ex: Tau = [1 0.5 0.33 0.25 0.2......1/n]. I can partially fix this issue by changing the initial values, but the parameters do not seem to be updating properly/organically.

How can I fix this issue?

Here is my code:

%Gen Data
N=25;
inv_ep = normrnd(.4,.1,[N,1]);
data.Ep1 = inv_ep(1:end-1,1);
.
.
.
opts = optimset ('Display','iter','TolX',2e-9,'TolFun',2e-9,'MaxIter',100000,'MaxFunEvals',10000000,'FunValCheck','off');
[param,fval,exitflag,output,lambda,grad,hessian] = fmincon(@(param) Posterior3(data,idx,param),param0,A,b,Aeq,beq,lb,ub,nonlcon,opts); 
err  = sqrt(diag(inv(hessian)));

function LL = Posterior3(data,idx,param)
%n = 13;  %number of obs: simple case
n = 24; %Will update to go by ID
% define components of the function:
Sigma_sq  = param(idx.Sigma_sq); % Variance of signal
ShockInfo = param(idx.ShockInfo); % Precision of the signal: n/Sigma_sq
Tau0      = param(idx.Tau);       % Variance of the prior mean
% Define vector to hold the epsilon shocks:
EpNew = zeros(n, 1);
% Define a vector to hold the precision of the prior:
Tau      = zeros(n+1,1);
Tau(1,:) = Tau0(1,:); % Variance of prior going into period 1
for i=1:n % Start the sampling
    EpNew(i,:) = data.Ep2(i,:) + (data.Ep2(i,:) - data.Ep1(i,:)).*((1./ShockInfo).*(1./Tau(i,:))); %Ep 1.5
    
    Tau(i+1,:) = ((1./ShockInfo).*Tau(i,:))./(Tau(i,:) + 1./ShockInfo);
end
%    Tau = Tau(1:end-1);
%    check = ((1./Tau(:))./(1./Tau(:) + ShockInfo(:))).*data.Ep1 + (ShockInfo(:)./(1./Tau(:) + ShockInfo(:))).*EpShock(:);
EpShock = data.Ep1 - EpNew;
LL = -sum(log(normpdf(EpShock, 0, Sigma_sq)));
end
```

Answer 1

[Matt J]

If the exitflag says a minimum was found then there is probably nothing suspect about the solution or the way it was reached. The likely culprits are either the initial guess param0, or the input problem data you're feeding to fmincon has a bug and doesn't specify the problem you actually intend.

One thing though is It's probably not a good idea to express the loglikelihood calculation this way:

log(normpdf(EpShock, 0, Sigma_sq))

The log of normpdf simplifies to an expression without exponentials. You should probably implement that expression directly, since exponential operations can be numerically tricky.

[Matt J]

Another possibility is your tolerance settings. We have no way of knowing whether 'TolX',2e-9,'TolFun',2e-9 are appropriate without being given the means to run the optimization ourselves.

[William Rose]

Two great suggesitons from @Matt J.

Answer 2

[William Rose]

@Joshua Scott,

Here is what I understand from reading your code. Correct me if I am wrong.

Function Posterior3(data, idx, param) returns the negative of the log likelihood. You want to minimize Posterior3, i.e. maximize likelihood, by adjusting param. The vector param contains SigmaSq, ShockInfo, and Tau0. Function Posterior3 uses SigmaSq "as is" to compute the likelihood. However, ShockInfo and Tau0 and not used directly in the likelihood calculation. Tau0 is the starting point for the calculation of vector Tau (25x1). ShockInfo is used in the calculation of vector EpNew (24x1). EpNew is subtracted from vector data.Ep1 to get vector EpShock (24x1). The likelihood is the joint probability of the 24 nromally distributed values EpSHock, each of which is normally distributed with mean 0 and variance SigmaSq (scalar or vector, I'm not sure which).

Yousaid "as Tau updates, its value equals some variant of 1/n. Ex: Tau = [1 0.5 0.33 0.25 0.2......1/n]. I can partially fix this issue by changing the initial values, but the parameters do not seem to be updating properly/organically." Tau is internal to Posterior3(). Tau is not an argument to Posterior3 (but Tau0 is), and therefore Tau is not (directly) optimized by fmincon(). The va;ues you specified for vector Tau are exactly what you expect if Tau0=1 and ShockInfo=1, because in that case, the formula for Tau(i+1) simplifies to

    Tau(i+1) = Tau(i)/(Tau(i) + 1);

Unrecognized function or variable 'Tau'.

which produces the sequence Tau=[1, .500, .333, .250, .200, .167,...] If that sequence suprises you, then reconsider the formula for Tau(i+1) in the for loop.

Thoughts:

• fmincon(), as you use it, adjusts Tau0, SigmaSq, and ShockInfo. fmincon() does not adjust the values of vector Tau - at least not directly. You can compute vector Tau after the fact, based on the optimum value of Tau0 returned by fmincon().
• Perhaps fmincon() is having trouble because it is working in a high dimensional search space. See quesiton below about the dimensions of SigmaSq, ShockInfo, and Tau0.
• Do not include values in vector param if they do not affect the value of Posterior3. Doing so will cause fmincon() to waste time searching fruitlessly for better values of parameters that make no difference. See question below abuout whther Tau0 is a vector.
• In your posting, you refer to the unexpected values of Tau. But Tau is not an argument to Posterior3. Did nyou mean to refer to Tau0?

Questions:

• Are SigmaSq, ShockInfo, and Tau0 scalars or vectors? If they are vectors, what is their size?
• If Tau0 is a vector, then note that only the first element of Tau0 is used in the calculation of LL.
• Do Tau0 and ShockInfo have reciprocal units, such as seconds and 1/seconds? I hope so, because of how they are combined in the code.
• Is it possible that the values of vector Tau and vector EpNew, obtained from the for loop, will tend to have similar final values, regardless of what "starting point" Tau0 is used? If so, then Tau0 is not a good candidate for optimization, since it has little effect on the final likelihood.

Suggestions:

Since Tau is a 25x1 vector, your code

Tau(1,:) = Tau0(1,:); % Variance of prior going into period 1

can be written as

Tau(1) = Tau0(1,:); % Variance of prior going into period 1

which is more clear. If Tau0 is a scalar (as I suspect), then you can simplify even more, as follows:

Tau(1) = Tau0; % Variance of prior going into period 1

Since EpNew is a 24x1 vector, and Tau is a 25x1 vector, your code

for i=1:n % Start the sampling
    EpNew(i,:) = data.Ep2(i,:) + (data.Ep2(i,:) - data.Ep1(i,:)).*((1./ShockInfo).*(1./Tau(i,:))); %Ep 1.5
    Tau(i+1,:) = ((1./ShockInfo).*Tau(i,:))./(Tau(i,:) + 1./ShockInfo);
end

can be wrtten more clearly as

for i=1:n % Start the sampling
    EpNew(i) = data.Ep2(i,:) + (data.Ep2(i,:) - data.Ep1(i,:)).*((1./ShockInfo)*(1/Tau(i))); %Ep 1.5
    Tau(i+1) = ((1./ShockInfo)*Tau(i))./(Tau(i) + 1./ShockInfo);
end

If Data.Ep1 and Data.Ep2 are column vectors, as I suspect, then you can simplify further, as shown below:

for i=1:n % Start the sampling
    EpNew(i) = data.Ep2(i) + (data.Ep2(i) - data.Ep1(i))*((1./ShockInfo)*(1/Tau(i))); %Ep 1.5
    Tau(i+1) = ((1./ShockInfo)*Tau(i))./(Tau(i) + 1./ShockInfo);
end

If ShockInfo is a scalar, then you can simplify further, as shown below:

for i=1:n % Start the sampling
    EpNew(i) = data.Ep2(i) + (data.Ep2(i) - data.Ep1(i))/(ShockInfo*Tau(i)); %Ep 1.5
    Tau(i+1) = 1/(ShockInfo + 1/Tau(i));
end

You said Tau follows the sequence 1,.5, .33, .25, ... That is exactly what you expect, based on the formula for Tau(i) inside Posterior3, if ShockInfo=1 and Tau0=1.