001 function [LL, cov] = mdist2dlike(params,data1,data2,dist)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033 if nargin < 3,
034 error('Requires at least FOUR input arguments');
035 end
036
037 if (nargin< 4)|isempty(dist),
038 error('Too few inputs')
039 else
040 HDIST=lower(dist{2});
041 VDIST=lower(dist{1});
042 end
043
044
045
046 switch VDIST(1:2),
047 case 'ra', nv=1;
048 otherwise, nv=2;
049 end
050 switch HDIST(1:2),
051 case 'ra', nh=1;
052 otherwise, nh=2;
053 end
054
055 if nv+nh+1~=length(params)
056 error('param is not the right size')
057 end
058
059 data1=data1(:);
060 data2=data2(:);
061 [n, m] = size(data1);
062 [n2, m2] = size(data2);
063 if n~=n2
064 error('data1 and data2 must have equal size')
065 end
066
067 if nargout == 2 & max(m,n) == 1
068 error('To compute the 2nd output, the 2nd input must have at least two elements.');
069 end
070 phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
071 phat.dist=dist;
072
073 x = mdist2dpdf(data1,data2,phat,0)+eps;
074 LL = -sum(log(x));
075
076 if nargout > 1
077 np=nv+nh+1;
078 sparam=params;
079 delta = eps^.4;
080 delta2=delta^2;
081
082 dist0 ='mdist2dpdf';
083
084
085
086
087 H = zeros(np);
088 for ix=1:np,
089 sparam = params;
090 sparam(ix)= params(ix)+delta;
091 phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
092 x = feval(dist0,data1,data2,phat)+eps;
093 fp = sum(log(x));
094 sparam(ix) = params(ix)-delta;
095 phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
096 x = feval(dist0,data1,data2,phat)+eps;
097 fm = sum(log(x));
098 H(ix,ix) = (fp+2*LL+fm)/delta2;
099 for iy=ix+1:np,
100 sparam(ix) = params(ix)+delta;
101 sparam(iy) = params(iy)+delta;
102 phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
103
104 x = feval(dist0,data1,data2,phat)+eps;
105 fpp = sum(log(x));
106 sparam(iy) = params(iy)-delta;
107 phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
108 x = feval(dist0,data1,data2,phat)+eps;
109 fpm = sum(log(x));
110 sparam(ix) = params(ix)-delta;
111 phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
112 x = feval(dist0,data1,data2,phat)+eps;
113 fmm = sum(log(x));
114 sparam(iy) = params(iy)+delta;
115 phat.x{1}=params(1:nv);phat.x{2}=params(nv+1:end-1);phat.x{3}=params(end);
116 x = feval(dist0,data1,data2,phat)+eps;
117 fmp = sum(log(x));
118 H(ix,iy) = (fpp-fmp-fpm+fmm)/(4*delta2);
119 H(iy,ix) = H(ix,iy);
120 end
121 end
122
123 cov = -H\eye(np);
124
125 end
126