octave_packages/vrml-1.0.13/best_dir_cov.m

   1 ## Copyright (C) 2002 Etienne Grossmann <etienne@egdn.net>
   2 ##
   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
   6 ## version.
   7 ##
   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
  11 ## details.
  12 ##
  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/>.
  15
  16 ##       [cv,wx] = best_dir_cov(x,a,sx,wd)
  17 ##
  18 ## x    D x P     :
  19 ## a    P x W     : Same as in best_dir, but sx is compulsory.
  20 ## sx   P x 1     :
  21 ##
  22 ## wd (W+D) x 1   : ML estimate of [w;d]
  23 ##
  24 ## cv (W+D)x(W+D) : Covariance of the ML estimator at [w;d]
  25 ##
  26 ## wx (W+D)x(P*D) : derivatives of ML estimate wrt to observations
  27 ##
  28
  29 ## Author:        Etienne Grossmann <etienne@egdn.net>
  30 ## Last modified: Setembro 2002
  31
  32 function  [cv,wx] = best_dir_cov(x,a,sx,wd)
  33
  34 [D,P] = size (x);
  35 W = columns (a);
  36 WD = prod (size (wd));
  37
  38                                 # Check dimensions etc
  39 if prod(size(sx)) != P
  40   error ("sx has %d != %d elements", prod (size (sx)), P);
  41 end
  42 if WD != W+D
  43   error ("wd has %d != %d elements", WD, W+D);
  44 end
  45 if rows (a) != P
  46   error ("a has %d != %d rows", rows (a), P);
  47 end
  48 if any (sx <= 0)
  49   error ("sx has some nonpositive elements");
  50 end
  51
  52 sx = sx(:) ;
  53 wd = wd(:) ;
  54
  55 w = wd(1:W);
  56 d = wd(W+1:WD);
  57
  58 isig = diag(1./sx) ;            # Inverse of covariance matrix.
  59
  60                                 # All derivatives are 1/2 of true value.
  61
  62 dsw = [zeros(W,1);d];           # Derivative of constraint |d|^2=1
  63
  64                                 # Inverse of Hessian with side blocks
  65 #keyboard
  66 if 0,                           # Readable code, bigger matrices
  67   d2ww = inv([ [-a';x]*isig*[-a,x'], dsw ; dsw' , 0 ]) ;
  68
  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 ]) ;
  72 end
  73 ## if any(abs(D2ww(:)-d2ww(:))>sqrt(eps)),
  74 ##  printf("Whoa!! %g",max(abs(D2ww(:)-d2ww(:)))) ;
  75 ## end
  76
  77                                 # 2nd Derivatives wrt.  wd  and  x
  78
  79 ## d2wx = zeros(WD+1,D*P);        # (padded with a row of zeros)
  80 d2wx = zeros(WD,D*P);
  81                                 # Easy : wrt.  w  and  x
  82 d2wx(1:W,:) = - kron(d',((1./sx)*ones(1,W))'.*a') ;
  83
  84 x = x'(:) ;
  85
  86 y = eye(D);                     # tmp
  87 tmp = zeros(D,D*P) ;
  88                                 # wrt. d  and  x
  89 for i=1:D,
  90
  91   ## d2wx(W+i,(i-1)*P+1:i*P) = \
  92   ##    2*x'*(kron(y(i,:))
  93   d2wx(W+i,:) =                        \
  94       2*x'*kron(y(i,:),kron(d,isig)) - \
  95       w'*a'*isig*kron(y(i,:),eye(P)) ;
  96 end
  97
  98 ## wx = d2ww*d2wx ;
  99
 100 wx = d2ww(1:WD,1:WD)*d2wx(1:WD,:) ;
 101 cv = ((wx.*kron(ones(WD,D),sx'))*wx') ;
 102
 103 ## cv = (wx*kron(eye(D),isig)*wx')(1:WD,1:WD) ;
 104 #  if 0,
 105 #    cv = (wx*kron(eye(D),diag(sx))*wx')(1:WD,1:WD) ;
 106 #  elseif 0
 107 #    cv = ((wx.*kron(ones(WD+1,D),sx'))*wx')(1:WD,1:WD) ;
 108 #  end
 109 #  if any(abs(cv2(:)-cv(:))>sqrt(eps)),
 110 #    printf("whoa!! b_d_cov (2) : %f\n",max(abs(cv2(:)-cv(:))));
 111 #    keyboard
 112 #  end
 113
 114
 115
 116
 117
 118
 119
 120 endfunction
 121