Re: OT: finding the weights used in weighted least squares regression

[Ben Abbott]

On Apr 27, 2011, at 1:19 AM, Kamaraju S Kusumanchi wrote:

> Ben Abbott wrote:
>> 
>> If I understand what you'd like to do, I think the solution is ...
>> 
>> W = (A*(A\b)-b) ./ (A*x-b)
>> 
> 
> No. This is not the solution. Consider the following example
> 
> octave:10> n=5, m=2, W = diag(rand(n,1)), A = rand(n,m), b = rand(n,1), x = 
> (W*A) \ (W*b), (A*(A\b)-b) ./ (A*x-b)
> n =  5
> m =  2
> W =
> 
>   0.63718   0.00000   0.00000   0.00000   0.00000
>   0.00000   0.75466   0.00000   0.00000   0.00000
>   0.00000   0.00000   0.69982   0.00000   0.00000
>   0.00000   0.00000   0.00000   0.08992   0.00000
>   0.00000   0.00000   0.00000   0.00000   0.94621
> 
> A =
> 
>   0.833095   0.568790
>   0.792483   0.536521
>   0.702069   0.588282
>   0.036704   0.529085
>   0.238762   0.427252
> 
> b =
> 
>   0.86356
>   0.77943
>   0.77652
>   0.98000
>   0.49344
> 
> x =
> 
>   0.32941
>   0.97441
> 
> ans =
> 
>     2.07510
>    -7.79993
>     2.55568
>     0.23533
>   112.87186
> 
> 
> Hence W is not same as (A*(A\b)-b) ./ (A*x-b) .
> 
> thanks
> -- 
> Kamaraju S Kusumanchi
> http://malayamaarutham.blogspot.com/

Ok. I see my mistake now. Essentially you want to solve for W where ...

        diag (W) * (A * x - b) = c

But you don't know "c". However, you do know that "c" are the least square 
errors.

I don't see a direct solution. Would an iterative solution work?

I don't know how stable it would be but, is something like below acceptable?

        err = @(u) x - (diag (u) * A) \ (diag (u) * b);
        W = fsolve (err, ones (n, 1))

I expect the options for fsolve will need to be added to get an accurate result.

Ben