Re: leasqr problem

[Jose Rodriguez]

> On Sat, Jan 24, 2009 at 1:12 PM, Jose Rodriguez
> <address@hidden> wrote:
> > I'm trying to solve a problem with leasqr, where the observed
> > value, instead of being a set of points y=f(x), is the result of a
> > definite integral. I found that I can use leasqr if I adjust the
> > dimensions of the function to fit and its expected values to the
> > dimension of the x range:
> >
> > x=(0:0.1:pi);
> >
> > P=@(x)  1/2 * ( 3*cos(x).^2-1 );
> >
> > F=inline( " ones ( size(x,1), 1 ) * trapz( x, P(x).*
> > exp(-a*P(x)).*sin(x) ) ", "x", "a" );
> >
> > y=0.5 * ones ( size(x,1), 1 );
> > pin=1;
> >
> > [f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=leasqr(x, y, pin, F);
> >
> >
> > If I run this I obtain the value 'a' for which the integral:
> >
> > trapz( x, P(x) .* exp(-a*P(x)) .* sin(x) )
> >
> > is 0.5.
> >
> > My problem is that I also want to include a normalisation
> > condition, like:
> >
> > 1 = trapz( x, exp(-a*P(x)).*sin(x) )
> >
> > Could somebody give me some pointers as to how to proceed? Or am I
> > completely out of track?

Jordi Gutiérrez Hermoso wrote:

> Hi, I think you're a little muddled about how you're approaching the
> problem here.

Most definitely so :)

> You are trying to solve the overdetermined system
>
>      F(a) = 0.5
>      G(a) = 1
>
> where F(a) and G(a) are your two integrals, you only have one
> observation point, the vector [0.5;1]. So the best you can do is try
> to find the a s.t. the norm of [F(a); G(a)] - [0.5; 1] is minimal.
>
> In this case it's not so hard to just search.

See below.


Jaroslav Hajek wrote:

> I'm not sure I understand what you're trying to do.
> The condition 1 == trapz( x, exp(-a*P(x)).*sin(x) ) is generally
> enough to determine a,
> because you have a single equation with a single variable.
> So is the other equation,
> trapz( x, P(x) .* exp(-a*P(x)) .* sin(x) ) == 0.5, but generally, you
> can't satisfy a system of two equations
> simultaneously if you only have one variable. Such a system may be
> solved in a least-squares sense, but
> is that what you need?

My original explanation wasn't very clear. I actually have 
two variables, only that I called them a=[a(1) a(2)] to 
comply with what leasqr expects; i.e. a function f(x,p) 
where p can be written like above to fit multiple variables. 
The equations look better like this:

P(x)=some function of x;

known_value=trapz(x, a*exp(b*P(x)).*P(x).*sin(x));

1=trapz(x, a*exp(b*P(x)).*sin(x));

Being 'a' and 'b' what I want to find.

It's true that for this example I could do with some 
plotting, but this is the simplest situation I can have. The 
general case involves n Pn(x) that go in the exponential, 
each of them with its coefficient to fit and its 
'known_value' arising from the integral, plus the 
normalisation condition. The number of variables to fit 
always equals the number of equations I have.

Lest squares is what other people have used to solve this 
very same problem, but I don't know how to put it in Octave.


Regards

Jose