* @return the matrix
*/
public Matrix mysut_wgn( EDBUnit f, Matrix augMu, Matrix augSigma, double Q, Map<String, Integer> mapName ){
//f.print("");
Matrix ret = new Matrix(2);
// check if any variable has 0 variance.
Matrix DP = augSigma.diag() ; //Sigma 1.25 0 0 -> DP 1.25
// 0 5 0 5
// 0 0 1 1
// if 0 is not in DP, oi=0, otherwise, oi is the index of variables with 0 variance.
//oi = find_index(0, DP); //If there are evidences, oi has their index
Matrix oi = DP.find_index( 0 );
// index of hidden variables.
//hi = find(DP) ; // find( ): find indices of nonzero elements, so hi = {1, 2, 3}
Matrix hi = DP.find();
// if one or more variables are observed, or almost sure with variance close to 0.
//if ~isempty(oi) //setting observation values
// Sigma = Sigma(hi,hi) ;
// obs_mu = mu(oi) ;
// mu = mu(hi,1) ;
// obs = 1 ; // observations found
//else
// obs = 0 ; // everything is normal, no observations found
Matrix obs_mu = null;
if( !oi.isempty() ) {
obs_mu = augMu.getMatrix(oi, 0);
augSigma = augSigma.getMatrix(hi, hi);
augMu = augMu.getMatrix(hi, 0);
}
// Sigma = Sigma(hi,hi) ;
// obs_mu = mu(oi) ;
// mu = mu(hi,1) ;
// obs = 1 ; // observations found
//else
// obs = 0 ; // everything is normal, no observations found
// SUT parameters setting:
//alpha=1 ; beta=2; kappa=0 ;
//L = length(mu) ; // number of dimensions of prior ( 3 = parent 1 + parent 2 + this node )
//lambda = alpha^2*(L + kappa) - L ;
//gamma = sqrt(L + lambda) ;
//S = gamma * Sql(Sigma)' ; // squared sigma
//Sigma 1.25 0 0 -> S 1.9365 0 0
// 0 5 0 0 3.873 0
// 0 0 1 0 0 1.7321
Double alpha=1.0, beta=2.0, kappa=0.0, lambda=0.0, gamma = 0.0;
int L = augMu.getRowDimension();
lambda = (alpha*alpha)*(L + kappa) - L ;
gamma = Math.sqrt(L + lambda) ;
Matrix S = augSigma.sqr().times(gamma);
// sigma-points generation, (2*L + 1) sigma points in total.
// nsp = 2*L +1 ; //2*3+1 = 7
int nsp = 2*L +1 ;
// sp(hi,1) = mu ; //sp = -15
// 65
// 0
Matrix sp = new Matrix(hi, augMu);
// for i=2:L+1{
// sp(hi,i) = mu + S(:,i-1) ; //sp = -15 -13.064
// 65 65
// 0 0
// sp(hi,i+L) = mu - S(:,i-1) ; //sp = -15 -13.064 0 0 -16.936
// } 65 65 0 0 65
// 0 0 0 0 0
// second step
// -15 -13.064 -15 0 -16.936 -15
// 65 65 68.873 0 65 61.127
// 0 0 0 0 0 0
// third step
// -15 -13.064 -15 -15 -16.936 -15 -15
// 65 65 68.873 65 65 61.127 65
// 0 0 0 1.7321 0 0 -1.7321
for( int i = 1; i< L+1; i++){
sp.addToColumn(hi, i, augMu.plus( S.getMatrixByRow(i-1)) );
sp.addToColumn(hi, i+L, augMu.minus( S.getMatrixByRow(i-1)) );
}
//add the observations into sigma-point.
// repmat( [3], 1, 2 ) => [3]
// [4] [4]
// [3]
// [4]
//if( obs )
// sp(oi,:) = repmat(obs_mu, 1, nsp) ; //????
if( !oi.isempty() ) {
sp.addToRow(oi, new Matrix( sp.getColumnDimension(), obs_mu ) );
}
// associated weights for sigma-points
//tf = lambda/(L+lambda) ;
//wm = tf ; // wm is weight for mean.
//th = 0.5/(L+lambda) ;
//wm(1,2:2*L+1) = repmat(th, 1, 2*L) ; //0 0.16667 0.16667 0.16667 0.16667 0.16667 0.16667
//wc = wm ; // wc is weight for cov.
//wc(1,1) = tf + (1 - alpha^2 + beta) ; //2 0.16667 0.16667 0.16667 0.16667 0.16667 0.16667
double tf = lambda/(L+lambda);
double th = 0.5/(L+lambda);
Matrix wm = new Matrix(1, nsp, th);
wm.set(0, 0, tf);
Matrix wc = new Matrix(wm);
wc.set(0, 0, tf + (1 - alpha*alpha + beta));
//compute
//for k=1:nsp //nsp = 7
// S2 = sp(:,k); //-15
// //65
// //0
// post_sp(:,k) = subs(func, func_var, S2) ; //func = -2*SerumCalcium+BrainTumor+noise
//func_var = 'SerumCalcium' 'BrainTumor' 'noise'
Matrix post_sp = new Matrix(1, nsp);
for( int k = 0; k < nsp; k++ ){
//f.print("");
post_sp.set(0, k, substitution(0, f, sp.getMatrixByRow(k), mapName));
}
//mean after function transformation.
//post_mu = post_sp * wm' ; // post_sp = 95 91.127 98.873 96.732 98.873 91.127 93.268
// post_mu = 95
Matrix post_mu = wm.times(post_sp);
//diff = post_sp - post_mu ; // post_sp - 95
//diff2 = diff .^ 2 ; // 0 15 15 3 15 15 3
Matrix diff = post_sp.minus(post_mu);
Matrix diff2 = diff.arrayTimes(diff);
//P = diff2 * wc' ; //11
Matrix P = diff2.times(wc);
//if nargin < 5, Q=0; //????
//post.mu = post_mu ; // 95
//post.Sigma = P + Q ; //11, white Gaussian noise is additive.
ret.set(0,0, post_mu.get());
ret.set(1,0, P.get() + Q);
return ret;
}