How to implement a summation function for numerical solving?

classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

How to implement a summation function for numerical solving?

andrewcd
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).
Reply | Threaded
Open this post in threaded view
|

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

Liam Groener-2

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

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 = (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
Reply | Threaded
Open this post in threaded view
|

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

andrewcd
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])
Reply | Threaded
Open this post in threaded view
|

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

andrewcd
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.
Reply | Threaded
Open this post in threaded view
|

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

Kimmo Luoma
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
Reply | Threaded
Open this post in threaded view
|

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

John W. Eaton
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
Reply | Threaded
Open this post in threaded view
|

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

andrewcd
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)<AT;
                EV(n) = c(n);
                endif
        end

L = length(EV);
EV = stp*L+xmin
Reply | Threaded
Open this post in threaded view
|

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

Liam Groener-2
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)<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.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
Reply | Threaded
Open this post in threaded view
|

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

Liam Groener-2
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 = (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)<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.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:
> <PastedGraphic-1.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.8491e-07  -1.9658e-06

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?


_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

andrewcd
Thanks again Liam.  

I am running Octave 3.2.3 in Ubuntu (Linux), on a 64-bit 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_64-pc-linux-gnu".

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 least-biased 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.