gennoiseds.m

递归贝叶斯估计的工具包

    if Gaussian_flag
        Arg.type = 'combo-gaussian';
        Arg.tag  = Noise.tag;
        Arg.dim  = Noise.dim;
        Arg.noiseSources = Noise.noiseSources;
        NoiseDS = gennoiseds(Arg);
    else
        % Restructure NoiseDS data structure
        NoiseDS.type = 'NoiseDS';
        NoiseDS.ns_type = Noise.type;
        NoiseDS.cov_type = Noise.cov_type;
        NoiseDS.tag = Noise.tag;
        NoiseDS.N = Noise.N;
        NoiseDS.idxArr = Noise.idxArr;
        NoiseDS.dim = Noise.dim;
        NoiseDS.mu = Noise.mu;
        NoiseDS.cov = Noise.cov;
        NoiseDS.noiseSources = Noise.noiseSources;
        NoiseDS.sample = @sample_combo;
        NoiseDS.update = @update_combo;
        NoiseDS.likelihood = @likelihood_combo;
    end


%===============================================================================================
%--- Gaussian Mixture Model
case 'gmm'

    if isfield(ArgDS,'M')
        Noise.M = ArgDS.M;
    else
        error(' [ gennoiseds::gmm ] Number of mixture components not specified.');
    end

    %-- assign noise source mean vector
    if ~isfield(ArgDS,'mu')
        Noise.mu = zeros(Noise.dim,Noise.M);    % default value
    else
        if (size(ArgDS.mu)==[Noise.dim Noise.M])
            Noise.mu = ArgDS.mu;
        else
            error(' [ gennoiseds::gmm ] Centroid mean dimension error.');
        end
    end

    %-- check for and assign cov_type
    if isfield(ArgDS,'cov_type')
        Noise.cov_type = ArgDS.cov_type;
    else
        warning(' [ gennoiseds::gmm ] Covariance type field .cov_type not assigned!. Assuming default value, ''full''');
        Noise.cov_type = 'full';             % default cov_type
    end

    %-- assign mixing weights
    if ~isfield(ArgDS,'weights')
        Noise.weights = (1/Noise.M)*ones(1,Noise.M);
    else
        if (length(ArgDS.weights)==Noise.M)
            Noise.weights = ArgDS.weights/(sum(ArgDS.weights));   % assign normalized weights
        else
            error(' [ gennoiseds::gmm ] Incorrect number of mixing weights (priors).');
        end
    end


    %-- assign rest of noise source structure
    switch (Noise.cov_type)

    %.............................................................................................
    case {'full','diag'}

        %-- assign noise source covariance structure
        if ~isfield(ArgDS,'cov')
            Noise.cov = repmat(eye(Noise.dim),[1 1 Noise.M]);    % default value
            warning(' [ gennoiseds::gmm ] Covariance field .cov not assigned!. Assuming default unity value.');
        elseif ((Noise.M == 1) & (size(ArgDS.cov) == [Noise.dim Noise.dim])) | ...
               ((Noise.M  > 1) & (size(ArgDS.cov) == [Noise.dim Noise.dim Noise.M]))
            Noise.cov = ArgDS.cov;
        else
            error([' [ gennoiseds::gmm::full ] Noise source covariance matrix has incorrect dimensions.']);
        end

    %.............................................................................................
    case {'sqrt','sqrt-diag'}

        %-- assign noise source covariance structure
        if ~isfield(ArgDS,'cov')
            Noise.cov = repmat(eye(Noise.dim),[1 1 Noise.M]);    % default value
            warning(' [ gennoiseds::gmm ] Covariance field .cov not assigned!. Assuming default unity value.');
        elseif ((Noise.M == 1) & (size(ArgDS.cov) == [Noise.dim Noise.dim])) | ...
               ((Noise.M  > 1) & (size(ArgDS.cov) == [Noise.dim Noise.dim Noise.M]))
            Noise.cov = ArgDS.cov;
        else
            error([' [ gennoiseds::gmm::sqrt ] Noise source covariance matrix has incorrect dimensions.']);
        end

    %.............................................................................................
    otherwise
        error(' [ gennoiseds::gmm ] Unkown noise source cov_type.');
    end


    % Restructure NoiseDS data structure
    NoiseDS.type = 'NoiseDS';
    NoiseDS.ns_type = Noise.type;
    NoiseDS.cov_type = Noise.cov_type;
    NoiseDS.tag = Noise.tag;
    NoiseDS.dim = Noise.dim;
    NoiseDS.M = Noise.M;
    NoiseDS.weights = Noise.weights;
    NoiseDS.mu = Noise.mu;
    NoiseDS.cov = Noise.cov;
    NoiseDS.sample = @gmmsample;
    NoiseDS.likelihood = @likelihood_gmm;




%===============================================================================================

otherwise
  error([' [ gennoiseds ] Noise type ' type ' not supported.']);
end

NoiseDS.adaptMethod = [];



%***********************************************************************************************
%***                                                                                         ***
%***                               SUB FUNCTION BLOCK                                        ***
%***                                                                                         ***
%***********************************************************************************************

%===============================================================================================
function NoiseDS = update_gaussian(NoiseDS)

% Updates a Gaussian noise source
%
% This function is only a placeholder here since Gaussian noise sources are completely updated, i.e.
% mean and covariance are set externally, i.e. there are no INTERNAL structure to update.


%===============================================================================================
function noise = sample_gaussian(NoiseDS, N)

% Generate N samples of a noise source specified by the NoiseDS data structure

    switch NoiseDS.cov_type

    case 'full'
        A = chol(NoiseDS.cov)';
    case 'diag'
        A = diag(sqrt(diag(NoiseDS.cov)));
    case 'sqrt'
        A = NoiseDS.cov;
    case 'sqrt-diag'
        A = NoiseDS.cov;
    otherwise
        error(' [ sample_gaussian ] Unknown cov_type.');
    end

    noise = A*randn(NoiseDS.dim,N) + cvecrep(NoiseDS.mu,N);


%===============================================================================================
function llh = likelihood_combo_gaussian(NoiseDS, noise, idxVec)

% Calculates the likelihood of sample 'noise', given the noise model NoiseDS. If the optional index
% vector 'idxVec' is specified, only those sub-noise sources are used. The 'noise' vector's dimension should
% concur with the implied total dimensionality of 'idxVec'

    if (nargin == 2)
      idxVec = 1:NoiseDS.N;
    end

    numNS = length(idxVec);

    [dim,nov] = size(noise);
    idxArr = NoiseDS.idxArr;        % copy beginning/ending index array
    llh = ones(1,nov);

    % ... for each noise source
    for j=1:numNS,

        idx1 = idxArr(idxVec(j),1);
        idx2 = idxArr(idxVec(j),2);

        idxRange = idx1:idx2;

        dim = idx2-idx1+1;

        switch NoiseDS.cov_type

        case 'full'
            D  = det(NoiseDS.cov(idxRange,idxRange));
            iP = inv(NoiseDS.cov(idxRange,idxRange));
        case 'diag'
            D = prod(diag(NoiseDS.cov(idxRange,idxRange)));
            iP = diag(1./diag(NoiseDS.cov(idxRange,idxRange)));
        case 'sqrt'
            D = det(NoiseDS.cov(idxRange,idxRange))^2;
            iS = inv(NoiseDS.cov(idxRange,idxRange));
            iP = iS'*iS;
        case 'sqrt-diag'
            D = prod(diag(NoiseDS.cov(idxRange,idxRange)))^2;
            iP = diag(1./(diag(NoiseDS.cov(idxRange,idxRange)).^2));
        otherwise
            error(' [ likelihood_gaussian ] Unkown cov_type.');
        end

        X = noise - cvecrep(NoiseDS.mu(idxRange),nov);
        q = 1/sqrt((2*pi)^dim * D);

        llh = llh .* (q * exp(-0.5*diag(X'*iP*X)'));

    end

    llh = llh + 1e-99; % needed to avoid 0 likelihood (cause ill conditioning)


%===============================================================================================
function llh = likelihood_combo(NoiseDS, noise, idxVec)

% Calculates the likelihood of sample 'noise', given the noise model NoiseDS.
% 'idxVec' is an optional index vector that can be used to indicate which of the N sub-noise sources should be used.
% to calculate the likelihood... this also requires 'noise' to have the same dimension of the relevant sub-noise source.

    if (nargin == 2)
      idxVec = 1:NoiseDS.N;
    end

    numNS = length(idxVec);

    [dim,nov] = size(noise);
    idxArr = NoiseDS.idxArr;        % copy beginning/ending index array
    llh = ones(1,nov);

    % ... for each noise source
    for j=1:numNS,

        idx1 = idxArr(idxVec(j),1);
        idx2 = idxArr(idxVec(j),2);
        subNoiseDS = NoiseDS.noiseSources{idxVec(j)};
        llh = llh .* feval(subNoiseDS.likelihood, subNoiseDS, noise(idx1:idx2,:));

    end

    llh = llh + 1e-99; % needed to avoid 0 likelihood (cause ill conditioning)


%===============================================================================================
function noise = sample_combo(NoiseDS, N)

%  Generate N samples of a noise source specified by the NoiseDS data structure

    noise=zeros(NoiseDS.dim,N);     % setup noise sample output buffer

    idxArr = NoiseDS.idxArr;        % copy beginning/ending index array

    % ... for each noise source
    for j=1:NoiseDS.N
        subNoiseDS = NoiseDS.noiseSources{j};
        noise(idxArr(j,1):idxArr(j,2),:) = feval(subNoiseDS.sample, subNoiseDS, N);
    end



%===============================================================================================
function noise = sample_combo_gaussian(NoiseDS, N)

%  Generate N samples of a noise source specified by the NoiseDS data structure

    noise=cvecrep(NoiseDS.mu,N);    % setup noise sample output buffer

    idxArr = NoiseDS.idxArr;        % copy beginning/ending index array

    num = NoiseDS.N;

    for j=1:num
        ind1 = idxArr(j,1);
        ind2 = idxArr(j,2);
        switch NoiseDS.cov_type
        case 'full'
            A = chol(NoiseDS.cov(ind1:ind2,ind1:ind2))';
        case 'diag'
            A = diag(sqrt(diag(NoiseDS.cov(ind1:ind2,ind1:ind2))));
        case 'sqrt'
            A = NoiseDS.cov(ind1:ind2,ind1:ind2);
        case 'sqrt-diag'
            A = NoiseDS.cov(ind1:ind2,ind1:ind2);
        otherwise
            error(' [ sample_gaussian ] Unkown cov_type.');
        end
        noise(ind1:ind2,:) = noise(ind1:ind2,:) + A*randn(ind2-ind1+1,N);
    end



%===============================================================================================
function NoiseDS = update_combo_gaussian(NoiseDS)

% Updates a 'combination Gaussian' noise source which has N Gaussian sub noise sources. The global mean and covariance
% is updated externally and then this function is called to update the internal sub-noise source structure.
%

idxArr = NoiseDS.idxArr;

for j=1:NoiseDS.N,
    ind1 = idxArr(j,1);
    ind2 = idxArr(j,2);
    idxRange = ind1:ind2;
    NoiseDS.noiseSources{j}.mu = NoiseDS.mu(idxRange,1);
    NoiseDS.noiseSources{j}.cov = NoiseDS.cov(idxRange,idxRange);
end

%===============================================================================================
function noise = sample_gamma(NoiseDS, N)

% Generate N samples of a noise source specified by the NoiseDS data structure

    alpha = NoiseDS.alpha;
    beta = NoiseDS.beta;

    if (alpha==1)
        noise = -log(1-rand(1,N))*beta;
        return
    end

    flag=0;

    if (alpha<1)
        flag=1;
        alpha=alpha+1;
    end

    gamma=alpha-1;
    eta=sqrt(2.0*alpha-1.0);
    c=.5-atan(gamma/eta)/pi;

    y(N)=0;

    for k=1:N,
        aux=-.5;
        while(aux<0)
            y(k)=-.5;
            while(y(k)<=0)
                u=rand(1,1);
                y(k) = gamma + eta * tan(pi*(u-c)+c-.5);
            end
            v=-log(rand(1,1));
            aux=v+log(1.0+((y(k)-gamma)/eta)^2)+gamma*log(y(k)/gamma)-y(k)+gamma;
        end
    end

    if (flag==1)
        noise = y.*beta.*(rand(1,N)).^(1.0/(alpha-1));
    else
        noise = y.*beta;
    end



%===============================================================================================
function llh = likelihood_gamma(NoiseDS, noise)

% Calculates the likelihood of sample 'noise', given the noise model NoiseDS.

    llh = noise.^(NoiseDS.alpha-1) .* exp((-1/NoiseDS.beta)*noise);

%===============================================================================================
function llh = likelihood_gmm(NoiseDS, noise)

% Calculates the likelihood of sample 'noise', given the noise model NoiseDS.

[prior,likelihood,evidence] = gmmprobability(NoiseDS, noise);

llh = evidence;