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 nonworking 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) = (nmu).^2 end for n = 1:lx p(n) = e.^((lam1*f1(n) + lam2*f2(n))) end mu = ones(1:lx)*mu fcns(1) = (f1mu)'*p fcns(2) = (f2variance)'*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 6sided 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). 
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 nonworking > 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) = (nmu).^2 > end > > for n = 1:lx > p(n) = e.^((lam1*f1(n) + lam2*f2(n))) > end > > mu = ones(1:lx)*mu > > fcns(1) = (f1mu)'*p > fcns(2) = (f2variance)'*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 6sided > 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/Howtoimplementasummationfunctionfornumericalsolvingtp3593289p3593289.html > Sent from the Octave  General mailing list archive at Nabble.com. > _______________________________________________ > Helpoctave mailing list > [hidden email] > https://mailman.cae.wisc.edu/listinfo/helpoctave I'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 = (f1mu).^2; p = exp(lam1*f1 lam2*f2); mu(1:lx) = mu; fcns = [(f1mu)'*p,(f2variance)'*p]; end _______________________________________________ Helpoctave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/helpoctave 
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 = (f1mu).^2; for n = 1:lx; p(n) = e.^((lam1*f1(n) + lam2*f2(n))); end mu(1:lx) = mu; fcns = [(f1mu)*p', (f2variance)*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([xmin1,xmax+1,0,1]) 
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. 
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, 20110614 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/Howtoimplementasummationfunctionfornumericalsolvingtp3593289p3596124.html > Sent from the Octave  General mailing list archive at Nabble.com. > _______________________________________________ > Helpoctave mailing list > [hidden email] > https://mailman.cae.wisc.edu/listinfo/helpoctave _______________________________________________ Helpoctave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/helpoctave 
Administrator

In reply to this post by andrewcd
On 14Jun2011, 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 _______________________________________________ Helpoctave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/helpoctave 
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 = (xmu).^2; p = exp(lam1*f1  lam2*f2) ; fcns = [(f1mu)*p', (f2variance)*p']; end result = fsolve(@eqns, guess) p = e.^(result(1)*x  result(2)*(xmu).^2); Z = sum(p); P=p./Z; plot( x, P); sum(P) axis([xmin1,xmax+1,0,max(P)*2]) c = cumsum(P); for n = 1:length(c); if c(n)<AT; EV(n) = c(n); endif end L = length(EV); EV = stp*L+xmin 
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 = (xmu).^2; > p = exp(lam1*f1  lam2*f2) ; > fcns = [(f1mu)*p', (f2variance)*p']; > end > > result = fsolve(@eqns, guess) > > p = e.^(result(1)*x  result(2)*(xmu).^2); > > Z = sum(p); > P=p./Z; > > plot( x, P); > sum(P) > axis([xmin1,xmax+1,0,max(P)*2]) > > c = cumsum(P); > for n = 1:length(c); > if c(n)<AT; > EV(n) = c(n); > endif > end > > L = length(EV); > EV = stp*L+xmin >  guess = 0.50000 0.00000 result = 8.0663e08 2.5000e01 FVEC = 7.5485e06 7.7026e05 INFO = 2 OUTPUT = iterations = 185 successful = 184 funcCount = 213 ans = 1 EV = 2  guess = 0 0 result = 6.5087e09 2.5000e01 FVEC = 1.3681e06 1.6720e05 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) _______________________________________________ Helpoctave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/helpoctave PastedGraphic1.pdf (25K) Download Attachment 
On Jun 15, 2011, at 10:46 PM, Liam Groener wrote:
> 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 = (xmu).^2; >> p = exp(lam1*f1  lam2*f2) ; >> fcns = [(f1mu)*p', (f2variance)*p']; >> end >> >> result = fsolve(@eqns, guess) >> >> p = e.^(result(1)*x  result(2)*(xmu).^2); >> >> Z = sum(p); >> P=p./Z; >> >> plot( x, P); >> sum(P) >> axis([xmin1,xmax+1,0,max(P)*2]) >> >> c = cumsum(P); >> for n = 1:length(c); >> if c(n)<AT; >> 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.0663e08 2.5000e01 > > FVEC = > 7.5485e06 7.7026e05 > > INFO = 2 > OUTPUT = > iterations = 185 > successful = 184 > funcCount = 213 > > ans = 1 > EV = 2 >  > guess = > 0 0 > > result = > 6.5087e09 2.5000e01 > > FVEC = > 1.3681e06 1.6720e05 > > INFO = 2 > OUTPUT = > iterations = 70 > successful = 69 > funcCount = 90 > > ans = 1.0000 > EV = 2 > > The central part of the plot is: > <PastedGraphic1.pdf> > 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) guess = 1 1 result = 9.2657 6.9418 FVEC = 8.8491e07 1.9658e06 INFO = 1 OUTPUT = iterations = 79 successful = 78 funcCount = 95 ans = 1.00000 EV = 1.3000 The curve obviously has a sharper higher peak. Though different, this solution seems to be mathematically better. (Note that INFO = 1 in this case vs 2 previously.) It looks live your equations have more than one solution? _______________________________________________ Helpoctave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/helpoctave 
Thanks again Liam.
I am running Octave 3.2.3 in Ubuntu (Linux), on a 64bit lenovo thinkpad x201. When QtOctave starts, it tells me: GNU Octave, version 3.2.3 Copyright (C) 2009 John W. Eaton and others. This is free software; see the source code for copying conditions. There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. For details, type `warranty'. Octave was configured for "x86_64pclinuxgnu". What do you mean by tolerances? Is this the degree of change that is needed between each iteration before the equation is considered "solved"? How does one change these? And is it a change to fsolve or to some other command? I am pretty sure (but not perfectly sure) that the equation should only have one solution. Basically what the code is doing is taking statistical moment constraints and fitting the leastbiased distribution that satisfies those constraints. Once I solve my current problem I want to add in higher moments like skewness and kurtosis. In fact, once I add in skewness (f3 = ((x  mu).^3)./(variance.^1.5);), and I manually iterate by taking the results and plugging them into the guesses, I can often eventually end up at a stable, if wrong, solution. The EV is a sort of hacked numerical integral from inf to .5  it should match the mean if my code is working properly. Not only does EV consistently not match the mean, but skewness is backwards, and the solutions take several manual iterations to become stable. I think that I have both code problems (Octave not iterating enough) as well as some fundamental problems (my approach to implementing MaxEnt has some flaw in it). Any thoughts on how I can solve the former? Thoughts on the later would be a welcome bonus. 
Free forum by Nabble  Edit this page 