function [ EstModelParams ] = estimateparametersfromsubpop( ModelParams, AlgorithmParams, ModelSpec, maxNStatesForEstimates, DataSpec )
    % This function takes a sample of estimates of various model
    % parameters, possibly on different state dimensions, and possibly
    % weighted, and returns the relevant Monte Carlo posterior mean point
    % estimates for those parameters relevant to a particular state
    % dimension.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % maxNStatesForEstimates - the subpopulation we are interested in
    % computing parameter estimates for.
    %
    % We return parameter estimates on the 'true' state space, i.e. if we
    % used the augmented states for parameter inference, the relevant
    % parameter estimates such as the generator matrix will be converted to
    % 'true' states.
    %
    % This function effectively replaces
    % estimateparametersfromallsubpops(), which makes estimates based on
    % weighted samples from all subpopulations. It makes more sense to use
    % one subpopulation only, because those of different state dimension
    % effectively represent different models.
    %
    % We should estimate the number of states before calling this function.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 23.10.2014 - created. Copied code from
    % estimateparametersfromallsubpops.m.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 4,
        error('Insufficient arguments!')
    end
    if ~isfield(AlgorithmParams, 'nBlocks'),
        error('nBlocks must be a field of AlgorithmParams!')
    end
    if ~isfield(AlgorithmParams, 'nParticlesEachBlock'),
        error('nParticlesEachBlock must be a field of AlgorithmParams!')
    end
    if ~isfield(AlgorithmParams, 'nTrueStatesEachBlock'),
        error('nTrueStatesEachBlock must be a field of AlgorithmParams!')
    end
    if ~isfield(AlgorithmParams, 'logParticleWeights'),
        error('logParticleWeights must be a field of AlgorithmParams!')
    end
    if ~isfield(ModelSpec, 'maxNStates'),
        error('maxNStates must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'useAugStateModel'),
        error('useAugStateModel must be a field of ModelSpec!')
    end
    if ModelSpec.useAugStateModel,
        if ~isfield(ModelSpec, 'augStateComponents'),
            error('augStateComponents must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'stateSpaceTransformation'),
            error('stateSpaceTransformation must be a field of ModelSpec!')
        end
    end
    if ~isfield(ModelSpec, 'dominatingRateScalingFactor'),
        error('dominatingRateScalingFactor must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'spikeTrainModelIndicator'),
        error('spikeTrainModelIndicator must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'positionModelIndicator'),
        error('positionModelIndicator must be a field of ModelSpec!')
    end
    if (rem(maxNStatesForEstimates, 1) ~= 0) || (maxNStatesForEstimates < 1),
        error('maxNStatesForEstimates be a positive integer!')
    end
    if maxNStatesForEstimates > max(ModelSpec.maxNStates),
        error('maxNStatesForEstimates must be less than the maximum number of states used for inference!')
    end
    if length(ModelParams) ~= AlgorithmParams.nBlocks,
        error('ModelParams must be a struct array of length AlgorithmParams.nBlocks!')
    end
    if ModelSpec.modelIndicator == 2,
        if ~isfield(ModelParams, 'generatorMat'),
            error('generatorMat must be a field of ModelParams!')
        end
        if ~isfield(ModelParams, 'jumpRate'),
            error('jumpRate must be a field of ModelParams!')
        end
    end
    if ~isfield(ModelParams, 'logInitialStateDist'),
        error('logInitialStateDist must be a field of ModelParams!')
    end
    if ModelSpec.spikeTrainModelIndicator > 0,
        if nargin < 5,
            error('Insufficient arguments!')
        end
        if ~isfield(DataSpec, 'nSpikeTrains'),
            error('nSpikeTrains must be a field of DataSpec!')
        end
        if ~isfield(ModelParams, 'logSpikeRate'),
            error('logSpikeRate must be a field of ModelParams!')
        end
    end
    if ModelSpec.positionModelIndicator > 0,
        if ~isfield(ModelParams, 'driftRate'),
            error('driftRate must be a field of ModelParams!')
        end
        if ~isfield(ModelParams, 'asympMeanPosInd'),
            error('asympMeanPosInd must be a field of ModelParams!')
        end
        if (ModelSpec.positionJitterModelIndicator == 0) || (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
            if ~isfield(ModelParams, 'jitterCov'),
                error('jitterCov must be a field of ModelParams!')
            end
        end
%         if ModelSpec.positionModelIndicator == 1,
%             if ~isfield(ModelParams, 'posObsVar'),
%                 error('posObsVar must be a field of ModelParams!')
%             end
%         end
    end

    % Which blocks represent the state dimension of interest?
    blocksThisSubpopBinInds = AlgorithmParams.nTrueStatesEachBlock == maxNStatesForEstimates;
    blocksThisSubpopInds = find(blocksThisSubpopBinInds);
    if ~any(blocksThisSubpopBinInds),
        error('There are no particles representing models with the requested state dimension!')
    end
    logSampleWeightSum = sumlogprobmat(cell2mat(AlgorithmParams.logParticleWeights(blocksThisSubpopBinInds)), 1);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    if ModelSpec.modelIndicator == 1,
        
        logTransitionMatBlockSum = -inf([maxNStatesForEstimates, maxNStatesForEstimates]);
        for iBlock = blocksThisSubpopInds',
            if ModelSpec.useAugStateModel,
                % Get into true state parameterisation.
                trueStateBinInds = ModelSpec.augStateComponents(1:ModelSpec.stateSpaceTransformation(maxNStatesForEstimates), 2) == maxNStatesForEstimates;
                trueStatelogTransitionMat = ModelParams(iBlock).logTransitionMat(trueStateBinInds, trueStateBinInds, :);

            else
                trueStatelogTransitionMat = ModelParams(iBlock).logTransitionMat;
            end
            % Weighted partial sum: multiply by weight, sum over particles and
            % add to running total.
            logTransitionMatPartialSum = sumlogprobmat(bsxfun(@plus, trueStatelogTransitionMat, permute(AlgorithmParams.logParticleWeights{iBlock}, [3, 2, 1])), 3);
            logTransitionMatBlockSum = sumlogprobmat(cat(3, logTransitionMatBlockSum, logTransitionMatPartialSum), 3);
        end
        % Divide by the sum of weights on each state.
        EstModelParams.logTransitionMat = logTransitionMatBlockSum - logSampleWeightSum;

    elseif ModelSpec.modelIndicator == 2,
        % Jump rate and generator matrix.
        % For the generator, we look at the jump rates and transition rates
        % separately then reassemble them later. This is so we can take logs
        % without worrying about sign, for log product computations.
        logJumpRateBlockSum = -inf([maxNStatesForEstimates, 1]);
        logTransitionRatesBlockSum = -inf([maxNStatesForEstimates, maxNStatesForEstimates]);
        for iBlock = blocksThisSubpopInds',
            if ModelSpec.useAugStateModel,
                % Get into true state parameterisation.
                trueStateBinInds = ModelSpec.augStateComponents(1:ModelSpec.stateSpaceTransformation(maxNStatesForEstimates), 2) == maxNStatesForEstimates;
                % Take log for efficient multiplication. Remember opposite
                % sign later.
                trueStateLogJumpRate = log(-ModelParams(iBlock).jumpRate(trueStateBinInds, :));
                trueStatelogTransitionRates = ModelParams(iBlock).generatorMat(trueStateBinInds, trueStateBinInds, :);

            else
                trueStateLogJumpRate = log(-ModelParams(iBlock).jumpRate);
                trueStatelogTransitionRates = ModelParams(iBlock).generatorMat;
            end
            % Turn self-transition rates to 0.
            indexingVec = (1:maxNStatesForEstimates)' + maxNStatesForEstimates * (0:(maxNStatesForEstimates - 1))';
            indexingVec = bsxfun(@plus, indexingVec, maxNStatesForEstimates ^ 2 * (0:(AlgorithmParams.nParticlesEachBlock(iBlock) - 1)));
            trueStatelogTransitionRates(indexingVec(:)) = 0;
            % Take log for efficient multiplication.
            trueStatelogTransitionRates = log(trueStatelogTransitionRates);
            % Weighted partial sum: multiply by weight, sum over particles and
            % add to running total.
            logJumpRatePartialSum = sumlogprobmat(bsxfun(@plus, trueStateLogJumpRate, permute(AlgorithmParams.logParticleWeights{iBlock}, [2, 1])), 2);
            logJumpRateBlockSum = sumlogprobmat(cat(2, logJumpRateBlockSum, logJumpRatePartialSum), 2);

            logTransitionRatesPartialSum = sumlogprobmat(bsxfun(@plus, trueStatelogTransitionRates, permute(AlgorithmParams.logParticleWeights{iBlock}, [3, 2, 1])), 3);
            logTransitionRatesBlockSum = sumlogprobmat(cat(3, logTransitionRatesBlockSum, logTransitionRatesPartialSum), 3);
        end
        % Divide by the sum of weights on each state.
        logJumpRateBlockSum = logJumpRateBlockSum - logSampleWeightSum;
        logTransitionRatesBlockSum = logTransitionRatesBlockSum - logSampleWeightSum;
        % Correct jump rates.
        EstModelParams.jumpRate = -exp(logJumpRateBlockSum);
        % Assemble generator.
        EstModelParams.generatorMat = exp(logTransitionRatesBlockSum);
        indexingVec = (1:maxNStatesForEstimates) + maxNStatesForEstimates * (0:(maxNStatesForEstimates - 1));
        EstModelParams.generatorMat(indexingVec) = EstModelParams.jumpRate;
        % Dominating rate.
        EstModelParams.dominatingRate = max(abs(EstModelParams.jumpRate), [], 1)' * ModelSpec.dominatingRateScalingFactor;
        % Compute log transition probabilities (conditional on jumps of
        % dominating process).
        EstModelParams.logTransitionMat = formdominatingprocesstransitionmat( EstModelParams.generatorMat, EstModelParams.dominatingRate );
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Initial state distribution.
    if ModelSpec.useAugStateModel,
        EstModelParams.logInitialStateDist = -inf([maxNStatesForEstimates, 1]);
        EstModelParams.logInitialStateDist(1) = 0;
    else
        logInitialStateDistBlockSum = -inf([maxNStatesForEstimates, 1]);
        for iBlock = blocksThisSubpopInds',
            trueStatelogInitialStateDist = ModelParams(iBlock).logInitialStateDist;
            % Weighted partial sum: multiply by weight, sum over particles
            % and add to running total.
            logInitialStateDistPartialSum = sumlogprobmat(bsxfun(@plus, trueStatelogInitialStateDist, permute(AlgorithmParams.logParticleWeights{iBlock}, [2, 1])), 2);
            logInitialStateDistBlockSum = sumlogprobmat(cat(2, logInitialStateDistBlockSum, logInitialStateDistPartialSum), 2);
        end
        % Divide by the total sum of weights.
        EstModelParams.logInitialStateDist = logInitialStateDistBlockSum - logSampleWeightSum;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Spike train model parameters.
    if ModelSpec.spikeTrainModelIndicator > 0,
        logSpikeRateBlockSum = -inf([maxNStatesForEstimates, DataSpec.nSpikeTrains]);
        for iBlock = blocksThisSubpopInds',
            % Weighted partial sum: multiply by weight, sum over particles
            % and add to running total.
            logSpikeRatePartialSum = sumlogprobmat(bsxfun(@plus, ModelParams(iBlock).logSpikeRate, permute(AlgorithmParams.logParticleWeights{iBlock}, [3, 2, 1])), 3);
            logSpikeRateBlockSum = sumlogprobmat(cat(3, logSpikeRateBlockSum, logSpikeRatePartialSum), 3);
        end
        % Divide by the sum of weights on each state.
        EstModelParams.logSpikeRate = bsxfun(@minus, logSpikeRateBlockSum, logSampleWeightSum);
    end
    if ModelSpec.modelIndicator == 1,
        if ModelSpec.spikeTrainModelIndicator == 1,
            EstModelParams.logNonspikeRate = log(1 - exp(EstModelParams.logSpikeRate));
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Position model parameters.
    if ModelSpec.positionModelIndicator > 0,
        
        % Drift rate.
        driftRateBlockSum = -inf([maxNStatesForEstimates, 1]);
        for iBlock = blocksThisSubpopInds',
            % Weighted partial sum: multiply by weight, sum over particles
            % and add to running total.
            driftRatePartialSum = sumlogprobmat(bsxfun(@plus, log(ModelParams(iBlock).driftRate), permute(AlgorithmParams.logParticleWeights{iBlock}, [2, 1])), 2);
            driftRateBlockSum = sumlogprobmat(cat(2, driftRateBlockSum, driftRatePartialSum), 2);
        end
        % Divide by the sum of weights on each state.
        EstModelParams.driftRate = exp(bsxfun(@minus, driftRateBlockSum, logSampleWeightSum));

%         % Instantaneous drift rate.
%         EstModelParams.instDriftRate = -log(EstModelParams.driftRate) ./ DataSpec.posObsInterval;

        % Asymptotic mean.
        % We are going to take the MODE of the sample for this parameter,
        % since the mean is not going to be a linear position ind.
        allWeights = exp(cell2mat(AlgorithmParams.logParticleWeights(blocksThisSubpopBinInds)));
        asympMeanPosInd = cat(2, ModelParams(blocksThisSubpopBinInds).asympMeanPosInd)';
        stateInds = repmat(1:size(asympMeanPosInd, 2), [length(allWeights), 1]);
        summingInds = [asympMeanPosInd(:), stateInds(:)];
        weightEachPosInd = accumarray(summingInds, repmat(allWeights, [size(asympMeanPosInd, 2), 1]), [DataSpec.nValidDiscretePositions, size(asympMeanPosInd, 2)]);
        [~, EstModelParams.asympMeanPosInd] = max(weightEachPosInd, [], 1);
        EstModelParams.asympMeanPosInd = EstModelParams.asympMeanPosInd';
%         asympMeanPosIndBlockSum = -inf([maxNStatesForEstimates, 1]);
%         for iBlock = blocksThisSubpopInds',
%             % Weighted partial sum: multiply by weight, sum over particles
%             % and add to running total.
%             asympMeanPosIndPartialSum = sumlogprobmat(bsxfun(@plus, ModelParams(iBlock).asympMeanPosInd, permute(AlgorithmParams.logParticleWeights{iBlock}, [2, 1])), 2);
%             asympMeanPosIndBlockSum = sumlogprobmat(cat(2, asympMeanPosIndBlockSum, asympMeanPosIndPartialSum), 2);
%         end
%         % Divide by the sum of weights on each state.
%         EstModelParams.asympMeanPosInd = bsxfun(@minus, asympMeanPosIndBlockSum, logSampleWeightSum);
        
        % Jitter covariance.
        if (ModelSpec.positionJitterModelIndicator == 0) || (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
            jitterCovBlockSum = -inf([2, maxNStatesForEstimates]);
            for iBlock = blocksThisSubpopInds',
                % Weighted partial sum: multiply by weight, sum over particles
                % and add to running total.
                jitterCovPartialSum = sumlogprobmat(bsxfun(@plus, log(ModelParams(iBlock).jitterCov), permute(AlgorithmParams.logParticleWeights{iBlock}, [3, 2, 1])), 3);
                jitterCovBlockSum = sumlogprobmat(cat(3, jitterCovBlockSum, jitterCovPartialSum), 3);
            end
            % Divide by the sum of weights on each state.
            EstModelParams.jitterCov = exp(bsxfun(@minus, jitterCovBlockSum, logSampleWeightSum));

            % Total covariance.
            EstModelParams.posCov = EstModelParams.jitterCov;
        elseif ModelSpec.positionJitterModelIndicator == 3,
            jitterCovBlockSum = zeros([2, 2, maxNStatesForEstimates]);
            for iBlock = blocksThisSubpopInds',
                % Weighted partial sum: multiply by weight, sum over particles
                % and add to running total.
                % Problem is, off diagonal can be negative.
                negInds = ModelParams(iBlock).jitterCov < 0;
                ModelParams(iBlock).jitterCov(negInds) = -ModelParams(iBlock).jitterCov(negInds);
                [jitterCovPartialSum, sumNegInds] = sumlogprobmat2(bsxfun(@plus, log(ModelParams(iBlock).jitterCov), permute(AlgorithmParams.logParticleWeights{iBlock}, [4, 2, 3, 1])), 4, negInds);
                jitterCovBlockSum = jitterCovBlockSum + exp(jitterCovPartialSum);
            end
            % Divide by the sum of weights on each state.
            negInds = jitterCovBlockSum < 0;
            jitterCovBlockSum(negInds) = -jitterCovBlockSum(negInds);
            EstModelParams.jitterCov = exp(bsxfun(@minus, log(jitterCovBlockSum), logSampleWeightSum));
            EstModelParams.jitterCov(negInds) = -EstModelParams.jitterCov(negInds);

            % Total covariance.
            EstModelParams.posCov = EstModelParams.jitterCov;
        end
        
        if ModelSpec.modelIndicator == 1,
            EstModelParams.stationaryPosCov = EstModelParams.posCov;
        elseif ModelSpec.modelIndicator == 2,
%             [ EstModelParams.instDriftRate, EstModelParams.instJitterCov, EstModelParams.instPosObsVar, EstModelParams.stationaryPosCov ] = computeinstantaneousparameters( EstModelParams.driftRate, EstModelParams.jitterCov, EstModelParams.posObsVar, ModelSpec.positionModelIndicator, DataSpec.posObsInterval, ModelSpec.randomWalkState );
            [ EstModelParams.instDriftRate, EstModelParams.instJitterCov, EstModelParams.stationaryPosCov ] = computeinstantaneousparameters( EstModelParams.driftRate, EstModelParams.jitterCov, ModelSpec.positionModelIndicator, DataSpec.posObsInterval, ModelSpec.randomWalkState );
        end
        
    end
end