

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/helpoctaveI'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:
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave

Administrator

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
>
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:
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


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)
Well, I just tried a guess of [1,1]. In this case, I did get a considerably different result:
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.

