Thread Subject:

summation

Hi,

I try to calculate expected backorder amount wih matlab, however I cannot calculate by using for loop for summation. I try to write:

sigma(p*(x-S)P(X=x)) and x starts from S and goes to infinity. P(X=x) is probability mass function( poisson distributin)I can do it with loop, however I need to do it without loop because I will use it in my objective function, I will use fminunc. I hope you will help me. Thank you in advance.

Kind regards,

Ozge

"Ozge Taskan" <lordgy@yahoo.com> wrote in message <hsg4kv$3vg$1@fred.mathworks.com>...
> Hi,
>
> I try to calculate expected backorder amount wih matlab, however I cannot calculate by using for loop for summation. I try to write:
>
> sigma(p*(x-S)P(X=x)) and x starts from S and goes to infinity. P(X=x) is probability mass function( poisson distributin)I can do it with loop, however I need to do it without loop because I will use it in my objective function, I will use fminunc. I hope you will help me. Thank you in advance.
>
> Kind regards,
>
> Ozge

  Does that "sigma(p*(x-S)P(X=x))" symbol mean that S is an integer and you are summing the expression p*(x-S)*P(X=x) through all integers x from S to infinity where P(X=x) means the poisson probability value

 lambda^x * exp(-lambda) / x-factorial ,

and what is the significance of the 'p*' there? Is p just a constant multiplicative factor? Exactly what is it you are optimizing with 'fminunc'? How high a value are you expecting to give S?

  If you are doing the actual summation up to a sufficiently large value of x to approximate infinity, there may be a way of replacing that infinite sum with a finite one without that approximation.

  Can you please explain things more carefully? I (and probably a number of others in this group) are not familiar with the concept of "expected backorder".

  One last comment. There is nothing in 'fminunc' that prevents your objective function from being computed in a loop. After all, even vectorized algorithms are actually doing their own looping internally. You just have to worry about total computing time.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hshkd1$kb1$1@fred.mathworks.com>...
> "Ozge Taskan" <lordgy@yahoo.com> wrote in message <hsg4kv$3vg$1@fred.mathworks.com>...
> > Hi,
> >
> > I try to calculate expected backorder amount wih matlab, however I cannot calculate by using for loop for summation. I try to write:
> >
> > sigma(p*(x-S)P(X=x)) and x starts from S and goes to infinity. P(X=x) is probability mass function( poisson distributin)I can do it with loop, however I need to do it without loop because I will use it in my objective function, I will use fminunc. I hope you will help me. Thank you in advance.
> >
> > Kind regards,
> >
> > Ozge
>
> Does that "sigma(p*(x-S)P(X=x))" symbol mean that S is an integer and you are summing the expression p*(x-S)*P(X=x) through all integers x from S to infinity where P(X=x) means the poisson probability value
>
> lambda^x * exp(-lambda) / x-factorial ,
>
> and what is the significance of the 'p*' there? Is p just a constant multiplicative factor? Exactly what is it you are optimizing with 'fminunc'? How high a value are you expecting to give S?
>
> If you are doing the actual summation up to a sufficiently large value of x to approximate infinity, there may be a way of replacing that infinite sum with a finite one without that approximation.
>
> Can you please explain things more carefully? I (and probably a number of others in this group) are not familiar with the concept of "expected backorder".
>
> One last comment. There is nothing in 'fminunc' that prevents your objective function from being computed in a loop. After all, even vectorized algorithms are actually doing their own looping internally. You just have to worry about total computing time.
>
> Roger Stafford

Thank you Roger ,
Let me say it more clearly,the situation of item i is sent with driver j is shown with binary variable y(i,j) (1 if it is sent via driver j, 0 otherwise), therefore p(i,j) is the penalty cost of sending item i with driver j if item becomes late. My function is sum over i, j ==> p(i,j)*y(i,j)*(sum over x[(x-S(i))*P(X=x))] and P(X=x) is (lambda(i,j)*y(i,j))^x * exp(-lambda(i,j)*y(i,j)) / x-factorial and x starts from s(i) and goes to infinity or as you say may be to 1000. This part (x-S(i))*P(X=x)) is supposed to find an expected number. It is just one part of my objective function, in general I try to optimize my s(i) and y(i,j). So I try to assign drivers to items and find s(i) s for each item i. And I could not write this function to fminunc as a part of objective function.
thank you very much in advance.

Ozge

"Ozge Taskan" <lordgy@yahoo.com> wrote in message <hsr76f$pi6$1@fred.mathworks.com>...
> Thank you Roger ,
> Let me say it more clearly,the situation of item i is sent with driver j is shown with binary variable y(i,j) (1 if it is sent via driver j, 0 otherwise), therefore p(i,j) is the penalty cost of sending item i with driver j if item becomes late. My function is sum over i, j ==> p(i,j)*y(i,j)*(sum over x[(x-S(i))*P(X=x))] and P(X=x) is (lambda(i,j)*y(i,j))^x * exp(-lambda(i,j)*y(i,j)) / x-factorial and x starts from s(i) and goes to infinity or as you say may be to 1000. This part (x-S(i))*P(X=x)) is supposed to find an expected number. It is just one part of my objective function, in general I try to optimize my s(i) and y(i,j). So I try to assign drivers to items and find s(i) s for each item i. And I could not write this function to fminunc as a part of objective function.
> thank you very much in advance.
>
> Ozge

  I am not certain this will help, but you can replace that infinite sum with a finite expression. I refer to what you call:

 (x-S(i))*P(X=x)) ,

summed from x = S(i) to x = infinity. If you have the statistics toolbox, there is a function in it called 'poisscdf' which is the cumulative poisson distribution. The above sum can be expressed in terms of this function as follows:

 (lambda-s)*(1-poisscdf(s-2,lambda)) + ...
 s*lambda^(s-1)*exp(-lambda)/factorial(s-1) .

  However this formula should be used with some care. If s is so large that poisscdf(s-2,lambda) is very close to 1, then when it is subtracted from 1 in the expression it loses relative accuracy, and it would have been better to do the sum directly since only a few more terms would have sufficed to approach infinity appropriately in such cases.

  As far as the outer sum with respect to i and j is concerned, all the above quantities in the finite expression above are capable of receiving vector inputs, so the outer sum could also be done without any for-loops.

  By the way, I don't know why you included the binary y(i,j) in the above poisson expressions along with lambda(i,j), since that whole inner sum is multiplied by y(i,j).

Roger Stafford

Could you please explain
> (lambda-s)*(1-poisscdf(s-2,lambda)) + ...
> s*lambda^(s-1)*exp(-lambda)/factorial(s-1) .
part, did you get this from x*P(X)-S(i)*P(X)?
In addition, did you mean to say that by defining x,P(x) and S as in vector form, can I solve this without any for loop?

Thank you very much in advance.

Regards,
Ozge

"Ozge Taskan" <lordgy@yahoo.com> wrote in message <hss5jk$o48$1@fred.mathworks.com>...
> Could you please explain
> > (lambda-s)*(1-poisscdf(s-2,lambda)) + ...
> > s*lambda^(s-1)*exp(-lambda)/factorial(s-1) .
> part, did you get this from x*P(X)-S(i)*P(X)?
> In addition, did you mean to say that by defining x,P(x) and S as in vector form, can I solve this without any for loop?
>
> Thank you very much in advance.
>
> Regards,
> Ozge

  I'll try to explain that identity, Ozge. First, let us agree on a temporary notation just for this occasion. For 'lambda' I'll write 'b'. I'll use n! for factorial n. By

 sum{n=0:inf}#(b^n/n!)

I'll mean the sum:

 (b^0/0! + b^1/1! + b^2/2! + b^3/3! + ...) (infinitely many)

  The reasoning goes like this. Temporarily drop exp(-b) out of things.

sum{n=s:inf}#((n-s)*b^n/n!) =
sum{n=s:inf}#(n*b^n/n!) - sum{n=s:inf}#(s*b^n/n!) =
b*sum{n=s:inf}#(b^(n-1)/(n-1)!) - s*sum{n=s:inf}#(b^n/n!) =
b*sum{n=s-1:inf}#(b^n/n!) - s*(sum{n=s-1:inf}#(b^n/n!)-b^(s-1)/(s-1)!) =
(b-s)*(sum{n=0:inf)#(b^n/n!)-sum(n=0:s-2}#(b^n/n!)) + s*b^(s-1)/(s-1)! =
(b-s)*exp(b)*(1-poisscdf(s-2,b)) + s*b^(s-1)/(s-1)!

Now multiply both first and last by exp(-b) and you have your identity. (You perhaps see why I chose to use an abbreviated notation.)

  The answer to your second question is yes. Everything in the derived expression can handle vector arguments, so you should be able to use the 'sum' operator to sum over i and j and never use any for-loops. I'm not sure it will make it much faster but you can compare it with using for-loops and see.

  Don't forget my caution about using s values that are too large for this method and losing accuracy. If only Mathworks had a "backwards" cumulative distribution for poisson distributions, it would help in this issue, but I haven't heard of such a thing.

Roger Stafford

Roger Stafford wrote:
> "Ozge Taskan" <lordgy@yahoo.com> wrote in message
> <hss5jk$o48$1@fred.mathworks.com>...
>> Could you please explain > (lambda-s)*(1-poisscdf(s-2,lambda)) + ...
>> > s*lambda^(s-1)*exp(-lambda)/factorial(s-1) .
>> part, did you get this from x*P(X)-S(i)*P(X)?
>> In addition, did you mean to say that by defining x,P(x) and S as in
>> vector form, can I solve this without any for loop?
>>
>> Thank you very much in advance.
>>
>> Regards,
>> Ozge
>
> I'll try to explain that identity, Ozge. First, let us agree on a
> temporary notation just for this occasion. For 'lambda' I'll write
> 'b'. I'll use n! for factorial n. By
>
> sum{n=0:inf}#(b^n/n!)
>
> I'll mean the sum:
>
> (b^0/0! + b^1/1! + b^2/2! + b^3/3! + ...) (infinitely many)
>
> The reasoning goes like this. Temporarily drop exp(-b) out of things.
>
> sum{n=s:inf}#((n-s)*b^n/n!) = sum{n=s:inf}#(n*b^n/n!) -
> sum{n=s:inf}#(s*b^n/n!) =
> b*sum{n=s:inf}#(b^(n-1)/(n-1)!) - s*sum{n=s:inf}#(b^n/n!) =
> b*sum{n=s-1:inf}#(b^n/n!) - s*(sum{n=s-1:inf}#(b^n/n!)-b^(s-1)/(s-1)!) =
> (b-s)*(sum{n=0:inf)#(b^n/n!)-sum(n=0:s-2}#(b^n/n!)) + s*b^(s-1)/(s-1)! =
> (b-s)*exp(b)*(1-poisscdf(s-2,b)) + s*b^(s-1)/(s-1)!
>
> Now multiply both first and last by exp(-b) and you have your identity.
> (You perhaps see why I chose to use an abbreviated notation.)
>
> The answer to your second question is yes. Everything in the derived
> expression can handle vector arguments, so you should be able to use the
> 'sum' operator to sum over i and j and never use any for-loops. I'm not
> sure it will make it much faster but you can compare it with using
> for-loops and see.
>
> Don't forget my caution about using s values that are too large for
> this method and losing accuracy. If only Mathworks had a "backwards"
> cumulative distribution for poisson distributions, it would help in this
> issue, but I haven't heard of such a thing.
>
> Roger Stafford

While MathWorks does not have a backwards cdf for Poisson, you can use
identity

1-poisscdf(x,lambda) = gammainc(lambda,x+1,'lower')

This improves accuracy. -Ilya

Ilya Narsky <inarsky@mathworks.com> wrote in message <hsslrb$qjh$1@fred.mathworks.com>...
> Roger Stafford wrote:
> > "Ozge Taskan" <lordgy@yahoo.com> wrote in message
> > <hss5jk$o48$1@fred.mathworks.com>...
> >> Could you please explain > (lambda-s)*(1-poisscdf(s-2,lambda)) + ...
> >> > s*lambda^(s-1)*exp(-lambda)/factorial(s-1) .
> >> part, did you get this from x*P(X)-S(i)*P(X)?
> >> In addition, did you mean to say that by defining x,P(x) and S as in
> >> vector form, can I solve this without any for loop?
> >>
> >> Thank you very much in advance.
> >>
> >> Regards,
> >> Ozge
> >
> > I'll try to explain that identity, Ozge. First, let us agree on a
> > temporary notation just for this occasion. For 'lambda' I'll write
> > 'b'. I'll use n! for factorial n. By
> >
> > sum{n=0:inf}#(b^n/n!)
> >
> > I'll mean the sum:
> >
> > (b^0/0! + b^1/1! + b^2/2! + b^3/3! + ...) (infinitely many)
> >
> > The reasoning goes like this. Temporarily drop exp(-b) out of things.
> >
> > sum{n=s:inf}#((n-s)*b^n/n!) = sum{n=s:inf}#(n*b^n/n!) -
> > sum{n=s:inf}#(s*b^n/n!) =
> > b*sum{n=s:inf}#(b^(n-1)/(n-1)!) - s*sum{n=s:inf}#(b^n/n!) =
> > b*sum{n=s-1:inf}#(b^n/n!) - s*(sum{n=s-1:inf}#(b^n/n!)-b^(s-1)/(s-1)!) =
> > (b-s)*(sum{n=0:inf)#(b^n/n!)-sum(n=0:s-2}#(b^n/n!)) + s*b^(s-1)/(s-1)! =
> > (b-s)*exp(b)*(1-poisscdf(s-2,b)) + s*b^(s-1)/(s-1)!
> >
> > Now multiply both first and last by exp(-b) and you have your identity.
> > (You perhaps see why I chose to use an abbreviated notation.)
> >
> > The answer to your second question is yes. Everything in the derived
> > expression can handle vector arguments, so you should be able to use the
> > 'sum' operator to sum over i and j and never use any for-loops. I'm not
> > sure it will make it much faster but you can compare it with using
> > for-loops and see.
> >
> > Don't forget my caution about using s values that are too large for
> > this method and losing accuracy. If only Mathworks had a "backwards"
> > cumulative distribution for poisson distributions, it would help in this
> > issue, but I haven't heard of such a thing.
> >
> > Roger Stafford
>
> While MathWorks does not have a backwards cdf for Poisson, you can use
> identity
>
> 1-poisscdf(x,lambda) = gammainc(lambda,x+1,'lower')
>
> This improves accuracy. -Ilya

  Yes Ilya, that is a valid identity (provided x is an integer of course), and using 'gammainc' should give more accurate results in that formula I gave than 'poisscdf' when the latter is near 1. Thanks for suggesting it.

  Ozge, the new, improved formula for your sum, thanks to Ilya, is:

 (lambda-s)*gammainc(lambda,s) + lambda^s*exp(-lambda)/factorial(s-1)

and you won't need the statistics toolbox for this. This is also capable of being vectorized.

Roger Stafford

Minimize (∑_(i=1) till n)[(∑_(j=1) till 4)〖(∑_(x=S_i) till ∞) p_ij * y_ij *(x-S_i )P(X=x) 〗+(∑_(i=1) till n)〖(∑_(x=0) till (S_i)〖h_i *(S_i-x) 〗 P(X=x) 〗+(∑_(i=1) till n)(∑_(j=1) till 4)〖t_ij *y_ij *m_i 〗]

this is my objective function and now with your answer it becomes
p.*y.*(lambda.*y-s)*gammainc((lambda.*y),s) + (lambda.*y)^s*exp(-(lambda.*y))/factorial(s-1)+h.lambda.(1- gammainc(lambda.*y,s,'lower'))+t.*y.*m

now my question is that y is a binary variable and a decision variable, and s is also my decision variable. I try to optimize this function however according to my research, matlab cannot solve nonlinear integer problems with constraints. Can it be possible to solve such a problem in matlab?

thank you in advance.

Ozge