1 ## Copyright (C) 2002 Etienne Grossmann <etienne@egdn.net>
3 ## This program is free software; you can redistribute it and/or modify it under
4 ## the terms of the GNU General Public License as published by the Free Software
5 ## Foundation; either version 3 of the License, or (at your option) any later
8 ## This program is distributed in the hope that it will be useful, but WITHOUT
9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 ## You should have received a copy of the GNU General Public License along with
14 ## this program; if not, see <http://www.gnu.org/licenses/>.
16 ## [cv,wx] = best_dir_cov(x,a,sx,wd)
19 ## a P x W : Same as in best_dir, but sx is compulsory.
22 ## wd (W+D) x 1 : ML estimate of [w;d]
24 ## cv (W+D)x(W+D) : Covariance of the ML estimator at [w;d]
26 ## wx (W+D)x(P*D) : derivatives of ML estimate wrt to observations
29 ## Author: Etienne Grossmann <etienne@egdn.net>
30 ## Last modified: Setembro 2002
32 function [cv,wx] = best_dir_cov(x,a,sx,wd)
36 WD = prod (size (wd));
38 # Check dimensions etc
39 if prod(size(sx)) != P
40 error ("sx has %d != %d elements", prod (size (sx)), P);
43 error ("wd has %d != %d elements", WD, W+D);
46 error ("a has %d != %d rows", rows (a), P);
49 error ("sx has some nonpositive elements");
58 isig = diag(1./sx) ; # Inverse of covariance matrix.
60 # All derivatives are 1/2 of true value.
62 dsw = [zeros(W,1);d]; # Derivative of constraint |d|^2=1
64 # Inverse of Hessian with side blocks
66 if 0, # Readable code, bigger matrices
67 d2ww = inv([ [-a';x]*isig*[-a,x'], dsw ; dsw' , 0 ]) ;
69 else # Unreadable, smaller matrices
70 ## tmp = (1./sx)*ones(1,WD);
71 d2ww = inv( [ ([-a,x'].*((1./sx)*ones(1,WD)))'*[-a,x'], dsw ; dsw', 0 ]) ;
73 ## if any(abs(D2ww(:)-d2ww(:))>sqrt(eps)),
74 ## printf("Whoa!! %g",max(abs(D2ww(:)-d2ww(:)))) ;
77 # 2nd Derivatives wrt. wd and x
79 ## d2wx = zeros(WD+1,D*P); # (padded with a row of zeros)
82 d2wx(1:W,:) = - kron(d',((1./sx)*ones(1,W))'.*a') ;
91 ## d2wx(W+i,(i-1)*P+1:i*P) = \
94 2*x'*kron(y(i,:),kron(d,isig)) - \
95 w'*a'*isig*kron(y(i,:),eye(P)) ;
100 wx = d2ww(1:WD,1:WD)*d2wx(1:WD,:) ;
101 cv = ((wx.*kron(ones(WD,D),sx'))*wx') ;
103 ## cv = (wx*kron(eye(D),isig)*wx')(1:WD,1:WD) ;
105 # cv = (wx*kron(eye(D),diag(sx))*wx')(1:WD,1:WD) ;
107 # cv = ((wx.*kron(ones(WD+1,D),sx'))*wx')(1:WD,1:WD) ;
109 # if any(abs(cv2(:)-cv(:))>sqrt(eps)),
110 # printf("whoa!! b_d_cov (2) : %f\n",max(abs(cv2(:)-cv(:))));