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

[Liam Groener]

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?