# How to implement a summation function for numerical solving? Classic List Threaded 10 messages Open this post in threaded view
|

## How to implement a summation function for numerical solving?

 Hi All,  First post here, and apologies if this topic is searchable -- I am a relative newcommer to octave, and don't know the appropriate terms to search for. My problem: I need to solve a nonlinear system of equations containing a summation term running the length of a vector with length (lx).  I have no idea how to implement this.  I can make it work by brute force -- manually typing in the expansion of the summation.  Of course, this is not flexible or especially useful as (lx) changes.   What I want to do is solve for (lam1) and (lam2).  Here is my non-working code: ------------------------------------------------------------------------------ function fcns = eqns(z) mu = 3.5 variance = 2.9 lx = 6 lam1=z(1) lam2=z(2) for n=1:lx         f1(n) = n         end for n=1:lx         f2(n) = (n-mu).^2         end         for n = 1:lx         p(n) = e.^(-(lam1*f1(n) + lam2*f2(n)))         end mu = ones(1:lx)*mu fcns(1) = (f1-mu)'*p fcns(2) = (f2-variance)'*p end -------------------------------------------------------------------- What I want to do, but what Octave doesn't seem to let me do, is to make a vector (p) that contains (lam1) and (lam2), and then solves for them from fcns(1) and fcns(2). When I try to run this code, I get the following error: -------------------------------------------------------------------- error: operator -: nonconformant arguments (op1 is 1x6, op2 is 1x2x3x4x5x6) error: called from: error:   eqns at line 29, column 9 error:   /usr/share/octave/3.2.3/m/optimization/fsolve.m at line 177, column 6 error:   /home/andrew/Documents/test.m at line 36, column 8--------------------------------------------------------------------- In addition, scrolling up, I see that Octave calculated (p) as (p = 1 1 1 1 1 1).  This is NOT the vector of unsolved functions that I want to later solve. There must be a way to do it, but it probably involves some programming that I haven't learned yet.  I've spent a long time trying to figure this out, but I'm not sure where to start, and it is difficult to know which terms to search for, as I don't know the jargon yet. Any help is greatly appreciated. (Context: I am trying to teach myself how to program max ent using a 6-sided die as an example.  I'll apply this to some other stuff later -- this is just for the purpose of learning how to get the code working).
Open this post in threaded view
|

## Re: How to implement a summation function for numerical solving?

 On Jun 12, 2011, at 11:47 PM, andrewcd wrote: > Hi All,  First post here, and apologies if this topic is searchable -- I am a > relative newcommer to octave, and don't know the appropriate terms to search > for. > > My problem: > I need to solve a nonlinear system of equations containing a summation term > running the length of a vector with length (lx).  I have no idea how to > implement this.  I can make it work by brute force -- manually typing in the > expansion of the summation.  Of course, this is not flexible or especially > useful as (lx) changes.   > > What I want to do is solve for (lam1) and (lam2).  Here is my non-working > code: > ------------------------------------------------------------------------------ > function fcns = eqns(z) > mu = 3.5 > variance = 2.9 > lx = 6 > > lam1=z(1) > lam2=z(2) > > for n=1:lx > f1(n) = n > end > > for n=1:lx > f2(n) = (n-mu).^2 > end > > for n = 1:lx > p(n) = e.^(-(lam1*f1(n) + lam2*f2(n))) > end > > mu = ones(1:lx)*mu > > fcns(1) = (f1-mu)'*p > fcns(2) = (f2-variance)'*p > > end > -------------------------------------------------------------------- > > What I want to do, but what Octave doesn't seem to let me do, is to make a > vector (p) that contains (lam1) and (lam2), and then solves for them from > fcns(1) and fcns(2). > > When I try to run this code, I get the following error: > -------------------------------------------------------------------- > error: operator -: nonconformant arguments (op1 is 1x6, op2 is 1x2x3x4x5x6) > error: called from: > error:   eqns at line 29, column 9 > error:   /usr/share/octave/3.2.3/m/optimization/fsolve.m at line 177, column > 6 > error:   /home/andrew/Documents/test.m at line 36, column 8 > --------------------------------------------------------------------- > > In addition, scrolling up, I see that Octave calculated (p) as (p = 1 1 1 1 > 1 1).  This is NOT the vector of unsolved functions that I want to later > solve. > > There must be a way to do it, but it probably involves some programming that > I haven't learned yet.  I've spent a long time trying to figure this out, > but I'm not sure where to start, and it is difficult to know which terms to > search for, as I don't know the jargon yet. > > Any help is greatly appreciated. > > (Context: I am trying to teach myself how to program max ent using a 6-sided > die as an example.  I'll apply this to some other stuff later -- this is > just for the purpose of learning how to get the code working). > > -- > View this message in context: http://octave.1599824.n4.nabble.com/How-to-implement-a-summation-function-for-numerical-solving-tp3593289p3593289.html> Sent from the Octave - General mailing list archive at Nabble.com. > _______________________________________________ > Help-octave mailing list > [hidden email] > https://mailman.cae.wisc.edu/listinfo/help-octaveI'm not completely clear in what you are trying to do. However, looking at your code; I expect this is what you wanted: function fcns = eqns(z) mu = 3.5; variance = 2.9; lx = 6; lam1=z(1); lam2=z(2); f1 = 1:lx; f2 = (f1-mu).^2; p = exp(-lam1*f1 -lam2*f2); mu(1:lx) = mu; fcns = [(f1-mu)'*p,(f2-variance)'*p]; end _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave
Open this post in threaded view
|

## Re: How to implement a summation function for numerical solving?

 Thanks Liam!  You helped me get it working.  I figured out that I had another couple of bugs, since fixed.  For posterity, here is the code to calculate a maximum entropy distribution given a mean and a variance.  No guarantees that it is the most efficient way to do this, but it will suit my purposes for now. Code: ------------------------------------------------------------------------------------------------------------------------------- clear all mu = 2 xmin = -5 xmax = 10 x = xmin:.2:xmax lx = length(x) stp = .2 function fcns = eqns(z) mu = 2; variance = 2; xmin = -5; xmax = 10; stp = .2 x = xmin:stp:xmax; lx = length(x); lam1=z(1); lam2=z(2); f1 = xmin:stp:xmax; f2 = (f1-mu).^2; for n = 1:lx;         p(n) = e.^(-(lam1*f1(n) + lam2*f2(n)));         end mu(1:lx) = mu; fcns = [(f1-mu)*p', (f2-variance)*p']; end guess = [0,0 ]; result = fsolve(@eqns, guess); p = zeros(1,lx); for n=1:lx;         p(n) = e.^(-((result(1)*x(n))+result(2)*(x(n)-mu).^2));         end         Z = sum(p); P=p./Z plot( x, P); sum(P) axis([xmin-1,xmax+1,0,1])
Open this post in threaded view
|

## Re: How to implement a summation function for numerical solving?

 One more question: You have probably noticed that my code defines variables twice -- once inside and once outside of the function.  This is a pain, obviously.  But variables defined outside of the function don't show up inside the function, and vice versa.  There must be a way around this.  Any hints?  Thanks.
Open this post in threaded view
|

## Re: How to implement a summation function for numerical solving?

 Hello, I have not read your complete code, but if you define some variable global then you could do something like global x; function r=fcn(y) global x; ...Do something endfunction This way x is also defined inside function. Kimmo On Tue, 2011-06-14 at 03:02 -0700, andrewcd wrote: > One more question: > > You have probably noticed that my code defines variables twice -- once > inside and once outside of the function.  This is a pain, obviously.  But > variables defined outside of the function don't show up inside the function, > and vice versa.  There must be a way around this.  Any hints?  Thanks. > > -- > View this message in context: http://octave.1599824.n4.nabble.com/How-to-implement-a-summation-function-for-numerical-solving-tp3593289p3596124.html> Sent from the Octave - General mailing list archive at Nabble.com. > _______________________________________________ > Help-octave mailing list > [hidden email] > https://mailman.cae.wisc.edu/listinfo/help-octave_______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave
Open this post in threaded view
|

## Re: How to implement a summation function for numerical solving?

 Administrator In reply to this post by andrewcd On 14-Jun-2011, andrewcd wrote: | One more question: | | You have probably noticed that my code defines variables twice -- once | inside and once outside of the function.  This is a pain, obviously.  But | variables defined outside of the function don't show up inside the function, | and vice versa.  There must be a way around this.  Any hints?  Thanks. Pass more than one input as needed.  For example,   mu = 2   xmin = -5   xmax = 10   x = xmin:.2:xmax   lx = length(x)   stp = .2   function result = eqns (z, mu, xmin, xmax, x, lx, stp)     ...   endfunction Then you can set up the call to your function with   mu = 2   xmin = -5   xmax = 10   x = xmin:.2:xmax   lx = length(x)   stp = .2   z = something;   result = eqns (z, mu, xmin, xmax, x, lx, stp) If you need to pass this function to fsolve, which expects a function of a single variable, then you can use an anonymous function:   mu = 2   xmin = -5   xmax = 10   x = xmin:.2:xmax   lx = length(x)   stp = .2   result = fsolve (@(z) eqns (z, mu, xmin, xmax, lx, stp), guess) All variables not listed in the argument list of the anonymous function are picked up from the surrounding context.  So in   @(z) eqns (z, mu, xmin, xmax, lx, stp) everything but Z is set from the variables that are in scope at the point where the anonymous function is created. Does that help? jwe _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave
Open this post in threaded view
|

## Re: How to implement a summation function for numerical solving?

 Thanks everybody for all of your help so far.  Now I have a new problem -- fsolve is highly sensitive to the initial "guess" vector.  I guess that it is not iterating enough, or something.  Or, since it is guessing at the value of an exponent, small changes by its standards lead to big changes in my results.  Any thoughts on how I can get around this?  Increase the number of iterations somehow?  Or is there a different way to solve simultaneous nonlinear equations? To get a flavor for my problem, try to run my code below, using guess = [.5,0] and guess = [0,0].  The code *should* eventually spit out a normal curve with a mean and variance of 2.  But the wrong guess gives the wrong curve. Thanks again, and thanks in advance. Code: ------------------------------------------------------------------------------------------------ clear all global mu; global x; global lx; global xmin; global variance; AT = .5; xmin = -10; xmax = 10; stp = .1; x = xmin:stp:xmax; lx = length(x); mu = 2; variance = 2; guess = [.5,0] function fcns = eqns(z); global x; global lx; global mu; global variance; lam1=z(1); lam2=z(2); f1 = x; f2 = (x-mu).^2; p = exp(-lam1*f1 - lam2*f2) ; fcns = [(f1-mu)*p', (f2-variance)*p']; end result = fsolve(@eqns, guess) p = e.^(-result(1)*x - result(2)*(x-mu).^2); Z = sum(p); P=p./Z; plot( x, P); sum(P) axis([xmin-1,xmax+1,0,max(P)*2]) c = cumsum(P); for n = 1:length(c);         if c(n)
Open this post in threaded view
|

## Re: How to implement a summation function for numerical solving?

 On Jun 15, 2011, at 7:21 AM, andrewcd wrote: > Thanks everybody for all of your help so far.  Now I have a new problem -- > fsolve is highly sensitive to the initial "guess" vector.  I guess that it > is not iterating enough, or something.  Or, since it is guessing at the > value of an exponent, small changes by its standards lead to big changes in > my results.  Any thoughts on how I can get around this?  Increase the number > of iterations somehow?  Or is there a different way to solve simultaneous > nonlinear equations? > > To get a flavor for my problem, try to run my code below, using guess = > [.5,0] and guess = [0,0].  The code *should* eventually spit out a normal > curve with a mean and variance of 2.  But the wrong guess gives the wrong > curve. > > Thanks again, and thanks in advance. > > Code: > ------------------------------------------------------------------------------------------------ > clear all > > global mu; > global x; > global lx; > global xmin; > global variance; > AT = .5; > xmin = -10; > xmax = 10; > stp = .1; > x = xmin:stp:xmax; > lx = length(x); > mu = 2; > variance = 2; > > guess = [.5,0] > > function fcns = eqns(z); > global x; > global lx; > global mu; > global variance; > lam1=z(1); > lam2=z(2); > f1 = x; > f2 = (x-mu).^2; > p = exp(-lam1*f1 - lam2*f2) ; > fcns = [(f1-mu)*p', (f2-variance)*p']; > end > > result = fsolve(@eqns, guess) > > p = e.^(-result(1)*x - result(2)*(x-mu).^2); > > Z = sum(p); > P=p./Z; > > plot( x, P); > sum(P) > axis([xmin-1,xmax+1,0,max(P)*2]) > > c = cumsum(P); > for n = 1:length(c); > if c(n) EV(n) = c(n); > endif > end > > L = length(EV); > EV = stp*L+xmin > Hi, Since no one else has answered your question, I thought I'd take a crack at it. But, I don't know what your problem actually is. I ran the two cases you suggested with the default tolerance. The results are: ------------------------------------------------------- guess =    0.50000   0.00000 result =    8.0663e-08   2.5000e-01 FVEC =   -7.5485e-06   7.7026e-05 INFO =  2 OUTPUT =     iterations =  185     successful =  184     funcCount =  213 ans =  1 EV =  2 ----------------------------- guess =    0   0 result =   -6.5087e-09   2.5000e-01 FVEC =   -1.3681e-06   1.6720e-05 INFO =  2 OUTPUT =     iterations =  70     successful =  69     funcCount =  90 ans =  1.0000 EV =  2 The central part of the plot is: Here the the [.5,0] guess is plotted as a solid red line while the [0,0] guess results are overlaid in green dots. All in all, the agreement seems within the default tolerances to me. (You could tighten them if you wish.)  In responding, please let us know what version of Octave you are using and what operating system you are using. (I'm using octave 3.4 on Mac OS X) _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave PastedGraphic-1.pdf (25K) Download Attachment