How to optimise this equation?

Hey,

I have the following equation and I need to optimize this to find the best value for n:

Furthermore, I have the following data:

Tau = 5.2
E(t) = [0 0.0199867 0.0999334 0.1598934 0.1998668 0.1598934 0.1199201 0.0799467 0.05996 0.0439707 0.02998 0.011992 0]';
t = [0 1 2 3 4 5 6 7 8 9 10 12 14]';

Hopefully can someone help me to find the best way for optimizing this function, I have looked for some methods but didn't really knew how to apply them.

Thanks in advance,

Kind regards,

Danny

[Michal]

What is the optimized function? What is condition for optimal n? What is the definition of best value?

[Danny Helwegen]

What is the optimized function?

I don't really know what you want to hear?

What is condition for optimal n?

n should be a positive number between 1 and ~ 20

What is the definition of best value?

The lowest n tha'ts possible between 1 and ~20

[Michal]

So, do you mean by Optimization the best fit of the function E(t) on defined points, where n is optimized variable (1 <= n <= 20)? If yes, see answers. If no, please clarify.

As you can see, the optimized function in your case is not E(t), as you say, but

fun = @(n) norm(E - E_t(t,n,tau));

Answer 1

[Star Strider]

Sinc you are using the factorial function, you are limited to integer values for ‘n’.

Try this:

objfcn = @(n,t,tau) (n.^n .* t.^(n-1)) ./ (factorial(n-1) .* tau.^n) .* exp(-n.*t/tau);

for k = 1:20
    fcn(k) = norm(E - objfcn(k,t,Tau));
end
[Out,n] = min(fcn)

and the value of ‘n’ that produces the lowest residual norm ia 4.

If you want to use the gamma function instead:

objfcn = @(n,t) (n.^n .* t.^(n-1)) ./ (gamma(n) .* Tau.^n) .* exp(-n.*t/Tau);

[B,rsdnrm] = fminsearch(@(b) norm(E - objfcn(b,t)), 10)

evaluating to:

B =
     4.2811

[Danny Helwegen]

Hey Star, thanks for the reply.

One question, where does the 10 after norm(E - objfcn(b,t)), stand for?

[Star Strider]

My pleasure.

That is the initial parameter estimate for ‘n’. The fminsearch function (and all other nonlinear optimization functions that I am aware of) have to have a starting estimate for the parameters (here one parameter) the are required to estimate. I chose 10 here because it is in the centre of the range you specify.

Answer 2

[John D'Errico]

Must n be an integer? After all, you have factorial(n-1) in there. If it must be an integer then just test a list of integer values for n. That will not be a long list, as you will have massive numerical problems if n is large. Pick the value of n that makes you happiest.

Hint: no, if n is an integer, you will never find an integer value for n that will come even remorely close to that data.

Instead, you very much need to use gamma(n) in there, instead of factorial(n-1). They are the same for integers, but the gamma function extends factorial nicely onto the real line.

Next, just to be safe, I used realpow in this. Although the .^ operator should be fine here.

Tau = 5.2
E = [0 0.0199867 0.0999334 0.1598934 0.1998668 0.1598934 0.1199201 0.0799467 0.05996 0.0439707 0.02998 0.011992 0]';
t = [0 1 2 3 4 5 6 7 8 9 10 12 14]';

E_t = @(t,n,tau) realpow(n,n).*realpow(t,n-1).*exp(-n*t/tau)./gamma(n)./realpow(tau,n);

Having done that, it is worth testing your fucntion to see if it has a chance to fit your data. So I played around with a few arbitrary values for n, and 4.25 seems to work resonably well, except that it does not reach the peak.

fplot(@(t) E_t(t,4.25,5.2),[0,14]),hold on,plot(t,E,'o')

As you see, it fits all ot the points pretty well, but it misses that peak. However, if you increase n to the point where is is large enough to hit the peak, then it looks like utter crap for the rest of the data.

fplot(@(t) E_t(t,6,5.2),[0,14]),hold on,plot(t,E,'o')

Could I have used an optimization to find the best value of n? Well, yes. That would be easy now. But the best value of n will be a compromise, likely close to 4.25, as I found initially.

fun = @(n) norm(E - E_t(t,n,tau));
[nopt,fval,exitflag] = fminsearch(fun,5)

nopt =
           4.2811279296875
fval =
        0.0301968802494492
exitflag =
     1
     
fplot(@(t) E_t(t,nopt,5.2),[0,14]),hold on,plot(t,E,'o')     

This points out several things. First, that I have a good eye for getting the best fit, even without a computer. But it also suggests that your data does not fit that model perfectly, missing that peak, even when the rest of the data fots nicely..