Related to a climate statistics post here.

function [tau,mu]=roman3(x,h,w);

% Finds tau and mu to minimize the sum of squares

% \sum_n \sum_m h(m,n)*(x(m,n)-tau(n)-mu(m))^2,

% subject to the constraint \sum_m w(m)*mu(m)=0. Here x and h are MxN,

% and w is Mx1. Normally, one would have non-negative h, positive w, and

% \sum_m w(m)=1. Upon return, tau will be 1xN, mu will be Mx1.

% The method uses a Lagrange multiplier (lambda) and Gaussian elimination.

% The matrix h must not contain an all-zero row or column.

[M,N]=size(x); % x and h have M rows, N columns

A=zeros(M+N+1,M+N+1); % A will contain the linear system matrix

A(1:N,1:N+M)=[diag(sum(h,1)) h']; % dL/dtau, N equations

A(N+1:N+M,:)=[h diag(sum(h,2)) -w/2]; % dL/dmu, M equations

A(end,N+1:N+M)=w'; % dL/dlambda, one equation (the constraint)

b=[sum(h.*x,1)'; sum(h.*x,2); 0]; % b is the "right hand side"

y=A\b; % solve Ay=b for y, using Gaussian elim.

tau=y(1:N)';

mu=y(N+1:N+M); % note the multiplier lambda is always zero

return

Ann says: Man, you beat me to it. That's just what I was thinking!

## 0 comments:

Post a Comment