0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 2 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 3 % Compute log-likelihood
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 4 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 5 % M Nofrarias 15-06-09
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 6 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 7 % $Id: loglikelihood_matrix.m,v 1.5 2011/11/16 08:52:50 nikos Exp $
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 8 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 10 function [loglk snr]= loglikelihood_matrix(xn,in,out,noise,model,params,inModel,outModel)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 11 % parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 12 fs = in(1).objs(1).fs;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 13
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 14 % num. experiments
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 15 nexp = numel(in);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 16 % num. transfer functions
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 17 nmod = numel(model(1).objs(:));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 18 noutChannels = sqrt(numel(noise(1).objs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 19
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 20 % loop over experiments
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 21 loglk = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 22 snr = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 23 for i = 1:nexp
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 24 if ((numel(in(1).objs) == 1) && numel(out(1).objs) == 1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 25
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 26 freqs = in(i).objs(1).data.getX;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 27 % evaluate models
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 28 if(isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 29 h11 = model(1).objs(1).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 30 elseif (~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 31 h11 = outModel(1,1).y * model(1).getObjectAtIndex(1,1).setParams(params,xn).double';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 32 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 33
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 34 % spectra to variance
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 35 % (N*fs/2)* this multiplication is done now in mcmc
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 36 C11 = noise(1).objs(1).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 37
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 38 % compute elements of inverse cross-spectrum matrix
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 39 InvS11 = 1./C11;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 40
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 41 % compute log-likelihood terms first, all at once does not cancel the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 42 % imag part when multiplying x.*conj(x)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 43 in1 = in(i).objs(1).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 44 out1 = out(i).objs(1).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 45
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 46 tmplt1 = h11.*in1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 47
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 48 v1v1 = conj(out1 - tmplt1).*(out1 - tmplt1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 49
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 50 %computing SNR
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 51 snrexp = utils.math.stnr(tmplt1,0,out1,0,InvS11,0,0,0);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 52
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 53 snr = snr + 20*log10(snrexp);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 54
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 55 log1exp = sum(InvS11.*v1v1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 56
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 57 loglk = loglk + log1exp;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 58
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 59 elseif ((numel(in(1).objs) == 2) && numel(out(1).objs) == 2)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 60 freqs = in(i).objs(1).data.getX;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 61 % loop over models
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 62
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 63 if(isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 64 for j = 1:nmod
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 65 % evaluate models
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 66 h(:,j) = model(i).objs(j).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 67 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 68 elseif (~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 69 h(:,1) = outModel(1,1).y * model(i).getObjectAtIndex(1,1).setParams(params,xn).double';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 70 h(:,2) = outModel(2,1).y * model(i).getObjectAtIndex(1,1).setParams(params,xn).double';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 71 h(:,3) = outModel(1,2).y * model(i).getObjectAtIndex(2,2).setParams(params,xn).double';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 72 h(:,4) = outModel(2,2).y * model(i).getObjectAtIndex(2,2).setParams(params,xn).double';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 73 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 74
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 75 for j = 1:noutChannels^2
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 76 % spectra to variance
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 77 % (N*fs/2)* this multiplication is done now in mcmc
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 78 C(:,j) = noise(i).objs(j).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 79 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 80
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 81 % compute elements of inverse cross-spectrum matrix
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 82 detm = (C(:,1).*C(:,4) - C(:,2).*C(:,3));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 83 InvS11 = C(:,4)./detm; %1 4
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 84 InvS22 = C(:,1)./detm; %4 1
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 85 InvS12 = C(:,2)./detm; %2 2
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 86 InvS21 = C(:,3)./detm; %3 3
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 87
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 88 % compute log-likelihood terms first, all at once does not cancel the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 89 % imag part when multiplying x.*conj(x)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 90 in1 = in(i).objs(1).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 91 in2 = in(i).objs(2).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 92 out1 = out(i).objs(1).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 93 out2 = out(i).objs(2).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 94
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 95 tmplt1 = h(:,1).*in1 + h(:,3).*in2;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 96 tmplt2 = h(:,2).*in1 + h(:,4).*in2;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 97
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 98 % matrix index convention: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 99 v1v1 = conj(out1 - tmplt1).*(out1 - tmplt1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 100 v2v2 = conj(out2 - tmplt2).*(out2 - tmplt2);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 101 v1v2 = conj(out1 - tmplt1).*(out2 - tmplt2);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 102 v2v1 = conj(out2 - tmplt2).*(out1 - tmplt1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 103
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 104 %computing SNR
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 105 snrexp = utils.math.stnr(tmplt1,tmplt2,out1,out2,InvS11,InvS22,InvS12,InvS21);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 106
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 107 snr = snr + 20*log10(snrexp);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 108
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 109 log1exp = sum(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 110
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 111 loglk = loglk + log1exp;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 112
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 113 elseif ((numel(in(1).objs) == 4) && numel(out(1).objs) == 3)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 114 % here we are implementing only the magnetic case
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 115 % We have 4 inputs (the 4 conformator waveforms of the magnetic
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 116 % analysis and
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 117 % 3 outputs (that correspond to the IFO.x12 and IFO.ETA1 and
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 118 % IFO.PHI1
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 119
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 120
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 121 for j = 1:noutChannels^2
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 122 % spectra to variance
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 123
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 124 % (N*fs/2)* this factor multiplication is done now in mcmc,
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 125 % before splitting
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 126 C(:,j) = noise(i).objs(j).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 127 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 128 if( isempty(inModel) && ~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 129
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 130 freqs = in(i).objs(1).data.getX;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 131
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 132 % faster this way
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 133 h(:,1) = outModel(1,1).y * model(i).getObjectAtIndex(1,1).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 134 h(:,2) = outModel(2,1).y * model(i).getObjectAtIndex(1,1).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 135 h(:,3) = outModel(3,1).y * model(i).getObjectAtIndex(1,1).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 136 h(:,4) = outModel(1,1).y * model(i).getObjectAtIndex(1,2).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 137 h(:,5) = outModel(2,1).y * model(i).getObjectAtIndex(1,2).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 138 h(:,6) = outModel(3,1).y * model(i).getObjectAtIndex(1,2).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 139 h(:,7) = outModel(1,2).y * model(i).getObjectAtIndex(2,3).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 140 h(:,8) = outModel(2,2).y * model(i).getObjectAtIndex(2,3).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 141 h(:,9) = outModel(3,2).y * model(i).getObjectAtIndex(2,3).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 142 h(:,10) = outModel(1,3).y * model(i).getObjectAtIndex(3,4).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 143 h(:,11) = outModel(2,3).y * model(i).getObjectAtIndex(3,4).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 144 h(:,12) = outModel(3,3).y * model(i).getObjectAtIndex(3,4).setParams(params,xn).double;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 145
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 146
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 147 % compute elements of inverse cross-spectrum matrix
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 148 detm = (C(:,1).*C(:,5).*C(:,9) + ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 149 C(:,2).*C(:,6).*C(:,7) + ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 150 C(:,3).*C(:,4).*C(:,8) -...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 151 C(:,7).*C(:,5).*C(:,3) -...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 152 C(:,8).*C(:,6).*C(:,1) -...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 153 C(:,9).*C(:,4).*C(:,2));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 154
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 155
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 156 InvS11 = (C(:,5).*C(:,9) - C(:,8).*C(:,6))./detm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 157 InvS12 = -(C(:,4).*C(:,9) - C(:,7).*C(:,6))./detm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 158 InvS13 = (C(:,4).*C(:,8) - C(:,7).*C(:,5))./detm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 159 InvS21 = -(C(:,2).*C(:,9) - C(:,8).*C(:,3))./detm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 160 InvS22 = (C(:,1).*C(:,9) - C(:,7).*C(:,3))./detm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 161 InvS23 = -(C(:,1).*C(:,8) - C(:,7).*C(:,2))./detm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 162 InvS31 = (C(:,2).*C(:,6) - C(:,5).*C(:,3))./detm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 163 InvS32 = -(C(:,1).*C(:,6) - C(:,4).*C(:,3))./detm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 164 InvS33 = (C(:,1).*C(:,5) - C(:,4).*C(:,2))./detm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 165
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 166 % compute log-likelihood terms first, all at once does not cancel the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 167 % imag part when multiplying x.*conj(x)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 168 for ll = 1:noutChannels
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 169 outV(:,ll) = out(i).objs(ll).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 170 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 171 for kk = 1:model(i).ncols
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 172 inV(:,kk) = in(i).objs(kk).data.getY;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 173 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 174
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 175 % faster this way
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 176 v(:,1) = outV(:,1) - h(:,1).*inV(:,1) - h(:,4).*inV(:,2) - h(:,7).*inV(:,3) - h(:,10).*inV(:,4);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 177 v(:,2) = outV(:,2) - h(:,2).*inV(:,1) - h(:,5).*inV(:,2) - h(:,8).*inV(:,3) - h(:,11).*inV(:,4);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 178 v(:,3) = outV(:,3) - h(:,3).*inV(:,1) - h(:,6).*inV(:,2) - h(:,9).*inV(:,3) - h(:,12).*inV(:,4);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 179
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 180 v1v1 = conj(v(:,1)).*v(:,1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 181 v1v2 = conj(v(:,1)).*v(:,2);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 182 v1v3 = conj(v(:,1)).*v(:,3);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 183 v2v1 = conj(v(:,2)).*v(:,1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 184 v2v2 = conj(v(:,2)).*v(:,2);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 185 v2v3 = conj(v(:,2)).*v(:,3);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 186 v3v1 = conj(v(:,3)).*v(:,1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 187 v3v2 = conj(v(:,3)).*v(:,2);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 188 v3v3 = conj(v(:,3)).*v(:,3);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 189
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 190 log1exp = sum(InvS11.*v1v1 +...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 191 InvS12.*v1v2 +...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 192 InvS13.*v1v3 +...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 193 InvS21.*v2v1 +...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 194 InvS22.*v2v2 +...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 195 InvS23.*v2v3 +...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 196 InvS31.*v3v1 +...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 197 InvS32.*v3v2 +...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 198 InvS33.*v3v3);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 199
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 200 loglk = loglk + log1exp;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 201
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 202
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 203 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 204 error('For the magnetic case, implement an outModel and leave your inModel blank')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 205 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 206
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 207 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 208 error('Implemented cases: 1 input / 1output, 2 inputs / 2outputs (TN3045 analysis), and 4 inputs / 3 outpus (magnetic complete analysis model. Other cases have not been implemented yet. Sorry for the inconvenience)');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 209 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 210
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 211 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 212 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 213