function PriorParams = initialisepriordistributions( ModelSpec, DataSpec )
    % Form the prior distribution parameter matrices.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ModelSpec - struct with the following fields:
    % - maxNStates - positive integer; the upper limit on the (true) state
    % dimension we will consider.
    % - estimateNStates - logical scalar; if true we set a uniform prior
    % over number of states; if false we set a prior with probabilty one on
    % maxNStates for number of states.
    % - spikeTrainModelIndicator, positionModelIndicator - if either of
    % these are 0, it indicates we do not wish to use the respective model
    % (spike trains or position).
    %
    % DataSpec - struct containing the following fields:
    % haveSpikeTrainsData, havePositionData - logicals.
    % nSpikeTrains - positive integer, required if haveSpikeTrainsData is
    % true.
    % nValidDiscretePositions - the number of discrete positions that define the maze
    % area, required if havePositionData is true.
    %
    % PriorParams is a struct with the following fields:
    % - stateDimensionDist - a column vector normalised distribution (not
    % in the log domain) over possible maximum (true) state dimensions. For
    % the uniform prior, these probabilities should all be equal.
    % - transitionRatesDirichletHyperparams
    % - jumpRateGammaHyperparams - last two are to do with state process.
    % - spikeRateGammaHyperparams - spike rates.
    % - driftRateUniformHyperparams
    % - asympMeanPriorDist
    % - obsErrorPrecisionGammaHyperparams
    % - jitterPrecisionGammaHyperparams
    % - jitterCovPriorDof
    % - jitterCovPriorScaleMat - last six are position model parameters.
    %
    % Note: the prior scale matrix of the position model can be interpreted
    % as the sum of squared deviations cross-product matrix (SSDCPM) for
    % the average field of coverage of any state's distribution over space.
    % We take it to be the 'typical' 'state' field's sum of squared
    % deviations. The jitterPrecisionGammaHyperparams have a similar
    % interpretation.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 01.10.2014 - copied from old version. Updated to new model, including
    % priors on all new parameters. Removed replication of hyperparameters
    % over states when these are all the same - but kept old version in
    % (commented out) in case needed.
    % 09.10.2014 - added initial state distribution prior.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 1,
        error('Insufficient arguments!')
    end
    if ~isfield(ModelSpec, 'maxNStates'),
        error('maxNStates must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'estimateNStates'),
        error('estimateNStates must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'useAugStateModel'),
        error('useAugStateModel 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 ~isfield(ModelSpec, 'initialStateModelIndicator'),
        error('initialStateModelIndicator must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'jumpRatePriorIndicator'),
        error('jumpRatePriorIndicator must be a field of ModelSpec!')
    end
    if ModelSpec.spikeTrainModelIndicator > 0,
%         if ~isfield(DataSpec, 'nSpikeTrains'),
%             error('nSpikeTrains must be a field of DataSpec!')
%         end
        if ~isfield(ModelSpec, 'spikeRatePriorIndicator'),
            error('spikeRatePriorIndicator must be a field of ModelSpec!')
        end
    end
    if (~ModelSpec.useAugStateModel) && (ModelSpec.initialStateModelIndicator == 3),
        if ~isfield(ModelSpec, 'initialStateDistPriorIndicator'),
            error('initialStateDistPriorIndicator must be a field of ModelSpec!')
        end
    end
    if ModelSpec.estimateNStates,
        if ~isfield(ModelSpec, 'stateDimPriorIndicator'),
            error('stateDimPriorIndicator must be a field of ModelSpec!')
        end
    end
    if ModelSpec.positionModelIndicator > 0,
        if nargin < 2,
            error('Insufficient arguments!')
        end
        if ~isfield(ModelSpec, 'positionObsErrorPriorIndicator'),
            error('positionObsErrorPriorIndicator must be a field of ModelSpec!')
        end
        if ~isfield(DataSpec, 'nValidDiscretePositions'),
            error('nValidDiscretePositions must be a field of DataSpec!')
        end
        if ~isfield(DataSpec, 'discreteYDim'),
            error('discreteYDim must be a field of DataSpec!')
        end
        if ~isfield(DataSpec, 'discreteXDim'),
            error('discreteXDim must be a field of DataSpec!')
        end
        if ~isfield(DataSpec, 'spatialBinWidth'),
            error('spatialBinWidth must be a field of DataSpec!')
        end
        if ~isfield(ModelSpec, 'positionJitterModelIndicator'),
            error('positionJitterModelIndicator must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'jitterPriorIndicator'),
            error('jitterPriorIndicator must be a field of ModelSpec!')
        end
        if (ModelSpec.positionObsErrorPriorIndicator == 0) || (ModelSpec.positionObsErrorPriorIndicator == 3),
            if ~isfield(ModelSpec, 'positionObsErrorTolerance'),
                error('positionObsErrorTolerance must be a field of ModelSpec!')
            end
        end
        if ModelSpec.modelIndicator == 2,
            if ~isfield(ModelSpec, 'driftRatePriorLowerBound'),
                error('driftRatePriorLowerBound must be a field of ModelSpec!')
            end
            if ~isfield(ModelSpec, 'driftRatePriorUpperBound'),
                error('driftRatePriorUpperBound must be a field of ModelSpec!')
            end
        end
        if ~isfield(ModelSpec, 'stateFieldSizeAsProportionOfShortestMazeDim'),
            error('stateFieldSizeAsProportionOfShortestMazeDim must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'nPriorObservationsForPositionModel'),
            error('nPriorObservationsForPositionModel must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'priorSpatialCorrelationCoef'),
            error('priorSpatialCorrelationCoef must be a field of ModelSpec!')
        end
    end

    % State space dimension prior.
    if ~ModelSpec.estimateNStates,
        PriorParams.stateDimensionDist = zeros([ModelSpec.maxNStates, 1]);
        PriorParams.stateDimensionDist(end) = 1;
    else
        if (ModelSpec.positionModelIndicator > 0) && (ModelSpec.randomWalkState > 0),
            nStatesForInference = ModelSpec.maxNStates - 1;

        else
            nStatesForInference = ModelSpec.maxNStates;
        end
        if ModelSpec.stateDimPriorIndicator == 1,
            PriorParams.stateDimensionDist = 1 / nStatesForInference * ones([nStatesForInference, 1]);
        elseif ModelSpec.stateDimPriorIndicator == 2,
            PriorParams.stateDimensionDist = geopdf((1:nStatesForInference)', 1 / nStatesForInference);
            PriorParams.stateDimensionDist = PriorParams.stateDimensionDist ./ sum(PriorParams.stateDimensionDist);
        elseif ModelSpec.stateDimPriorIndicator == 3,
            PriorParams.stateDimensionDist = geopdf((1:nStatesForInference)', 1 / 2 / nStatesForInference);
            PriorParams.stateDimensionDist = PriorParams.stateDimensionDist ./ sum(PriorParams.stateDimensionDist);
        elseif ModelSpec.stateDimPriorIndicator == 4,
            PriorParams.stateDimensionDist = geopdf((1:nStatesForInference)', 1 / 4 / nStatesForInference);
            PriorParams.stateDimensionDist = PriorParams.stateDimensionDist ./ sum(PriorParams.stateDimensionDist);
        end
        if ModelSpec.randomWalkState > 0,
            PriorParams.stateDimensionDist = [0; PriorParams.stateDimensionDist];
        end
    end
    
    % Initial state distribution (only if we are not using augmented states
    % and actually wish to make inference for it).
    if (~ModelSpec.useAugStateModel) && (ModelSpec.initialStateModelIndicator == 3),
        if ModelSpec.initialStateDistPriorIndicator == 1
            PriorParams.initialStateDistDirichletHyperparams = ones([1, ModelSpec.maxNStates]);
        else
            error('Invalid initialStateDistPriorIndicator!')
        end
    end
    
%     PriorParams.transitionRatesDirichletHyperparams = ones([ModelSpec.maxNStates, ModelSpec.maxNStates - 1]);
    if ModelSpec.modelIndicator == 1,
        PriorParams.transitionRatesDirichletHyperparams = ones([1, ModelSpec.maxNStates]);
    elseif ModelSpec.modelIndicator == 2,
        % Generator matrix prior parameters.
        % For each of N state there are N - 1 transition rates required to
        % sample in prior. These have a Dirichlet prior, conditional on the
        % corresponding jump rate.
        PriorParams.transitionRatesDirichletHyperparams = ones([1, ModelSpec.maxNStates - 1]);

        % Jump rates have a Gamma prior.
        if ModelSpec.jumpRatePriorIndicator == 1,
    %         PriorParams.jumpRateGammaHyperparams = ones([ModelSpec.maxNStates, 2]);
            PriorParams.jumpRateGammaHyperparams = [1, 1];
        elseif ModelSpec.jumpRatePriorIndicator == 2,
            PriorParams.jumpRateGammaHyperparams = [1/2, 1/2];
        elseif ModelSpec.jumpRatePriorIndicator == 3,
            % Near-Jeffreys.
    %         PriorParams.jumpRateGammaHyperparams = repmat([1/2, 0.01], [ModelSpec.maxNStates, 1]);
            PriorParams.jumpRateGammaHyperparams = [1/2, 0.001];
        elseif ModelSpec.jumpRatePriorIndicator == 4,
            PriorParams.jumpRateGammaHyperparams = ModelSpec.jumpRateGammaHyperparams;
        else
            error('Invalid jumpRateGammaHyperparams parameter!')
        end
    end
    
    % Spike trains model spike rate prior parameters.
    if (ModelSpec.spikeTrainModelIndicator == 1) || (ModelSpec.spikeTrainModelIndicator == 2),
        if ModelSpec.spikeRatePriorIndicator == 1,
            % Gam(1, 1).
%             PriorParams.spikeRateGammaHyperparams = ones([ModelSpec.maxNStates, DataSpec.nSpikeTrains, 2]);
            PriorParams.spikeRateGammaHyperparams = [1, 1];
        elseif ModelSpec.spikeRatePriorIndicator == 2,
            % Gam(1/2, 1/2).
            PriorParams.spikeRateGammaHyperparams = [1/2, 1/2];
        elseif ModelSpec.spikeRatePriorIndicator == 3,
            % Near-Jeffreys prior.
%             hyperparams = [1/2, 0.01];
%             PriorParams.spikeRateGammaHyperparams = repmat(permute(hyperparams, [1, 3, 2]), [ModelSpec.maxNStates, DataSpec.nSpikeTrains, 1]);
            PriorParams.spikeRateGammaHyperparams = [1/2, 0.001];
        elseif ModelSpec.spikeRatePriorIndicator == 4,
            PriorParams.spikeRateGammaHyperparams = ModelSpec.spikeRateGammaHyperparams;
        elseif ModelSpec.spikeRatePriorIndicator == 5,
            % TEST.
%             hyperparams = [1/4, 1/2];
%             hyperparams = [5.7182e-001,  2.1432e-001];
%             PriorParams.spikeRateGammaHyperparams = repmat(permute(hyperparams, [1, 3, 2]), [ModelSpec.maxNStates, DataSpec.nSpikeTrains, 1]);
            PriorParams.spikeRateGammaHyperparams = [5.7182e-001,  2.1432e-001];
        else
            error('Invalid spikeRatePriorIndicator parameter!')
        end
    elseif ModelSpec.spikeTrainModelIndicator > 2,
        error('Invalid spike trains model parameter!')
    end

    % Position model prior parameters.
    if ModelSpec.positionModelIndicator == 1,
        % Drift rates.
        % Uniform lower and upper bounds.
        if ModelSpec.modelIndicator == 2,
            PriorParams.driftRateUniformHyperparams = [ModelSpec.driftRatePriorLowerBound, ModelSpec.driftRatePriorUpperBound];
        end
        
        % Asymptotic mean positions.
        % This is uniform over the discrete regions that define the maze.
        % We only use this for sampling from the prior, not the full
        % conditional, which is handled differently.
        PriorParams.asympMeanPriorDist = 1 / DataSpec.nValidDiscretePositions * ones([DataSpec.nValidDiscretePositions, 1]);
        
    elseif ModelSpec.positionModelIndicator == 2,
        % Drift rates.
        % Uniform lower and upper bounds.
        PriorParams.driftRateUniformHyperparams = [ModelSpec.driftRatePriorLowerBound, ModelSpec.driftRatePriorUpperBound];
        % For asymptotic mean position, use uniform distribution over
        % discrete positions.
        PriorParams.asympMeanPriorDist = 1 / DataSpec.nValidDiscretePositions * ones([DataSpec.nValidDiscretePositions, 1]);
    end
    % Position model 'jitter' model parameters (covariance matrix).
    if ModelSpec.positionModelIndicator > 0,
        if (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
            % Gamma prior for single dimension precision parameters.
            if ModelSpec.jitterPriorIndicator == 1,
                % Gam(1, 1).
                PriorParams.jitterPrecisionGammaHyperparams = [1, 1];
            elseif ModelSpec.jitterPriorIndicator == 2,
                % Gam(1/2, 1/2).
                PriorParams.jitterPrecisionGammaHyperparams = [1/2, 1/2];
            elseif ModelSpec.jitterPriorIndicator == 3,
                % Near-Jeffreys prior.
                PriorParams.jitterPrecisionGammaHyperparams = [1/2, 0.01];
            elseif ModelSpec.jitterPriorIndicator == 4,
                % Prior based on 'state field' ideal.
%                 stateFieldSsdev = ModelSpec.nPriorObservationsForPositionModel * min([DataSpec.discreteYDim, DataSpec.discreteXDim]) * ModelSpec.stateFieldSizeAsProportionOfShortestMazeDim * DataSpec.spatialBinWidth;
%                 stateFieldSsdev = ModelSpec.nPriorObservationsForPositionModel * (min([DataSpec.discreteYDim, DataSpec.discreteXDim]) * ModelSpec.stateFieldSizeAsProportionOfShortestMazeDim * DataSpec.spatialBinWidth) ^ 2;
                stateFieldSsdev = ModelSpec.nPriorObservationsForPositionModel * (min([DataSpec.discreteYDim, DataSpec.discreteXDim]) * ModelSpec.stateFieldSizeAsProportionOfShortestMazeDim * DataSpec.spatialBinWidth * DataSpec.posObsInterval) ^ 2;
%                 hyperparams = [ModelSpec.nPriorObservationsForPositionModel / 2, 2 / stateFieldSsdev];
%                 % Rate parameter of gamma for precision should be
%                 % reciprocal, as here, of what it would be for variance.
%                 % (Note reciprocal will be taken again for MATLAB's gamrnd
%                 % function which uses SCALE, not rate).
                hyperparams = [ModelSpec.nPriorObservationsForPositionModel / 2, stateFieldSsdev / 2];
                PriorParams.jitterPrecisionGammaHyperparams = hyperparams;
            elseif ModelSpec.jitterPriorIndicator == 5,
                % Specified.
                PriorParams.jitterPrecisionGammaHyperparams = ModelSpec.jitterGammaHyperparams;
            else
                error('Invalid jitterPriorIndicator parameter!')
            end
        elseif (ModelSpec.positionJitterModelIndicator == 3) && (ModelSpec.positionModelIndicator == 2),
            % Inverse-Wishart prior parameters: degrees of freedom and
            % scale matrix. Based on 'state field' ideal.
            % This must be at least 2.
            PriorParams.jitterCovPriorDof = ModelSpec.nPriorObservationsForPositionModel;
            PriorParams.jitterCovPriorScaleMat = [
                1, 0;
                0, 1
            ];
%             stateFieldSsdev = PriorParams.jitterCovPriorDof * min([DataSpec.discreteYDim, DataSpec.discreteXDim]) * ModelSpec.stateFieldSizeAsProportionOfShortestMazeDim;
            % Jitter is over pixels, still, not discrete positions.
            stateFieldSsdev = PriorParams.jitterCovPriorDof * min([DataSpec.discreteYDim, DataSpec.discreteXDim]) * ModelSpec.stateFieldSizeAsProportionOfShortestMazeDim * DataSpec.spatialBinWidth;
            PriorParams.jitterCovPriorScaleMat = PriorParams.jitterCovPriorScaleMat .* stateFieldSsdev;
            PriorParams.jitterCovPriorScaleMat(1, 2) = sqrt(PriorParams.jitterCovPriorScaleMat(1, 1)) * sqrt(PriorParams.jitterCovPriorScaleMat(2, 2)) * ModelSpec.priorSpatialCorrelationCoef;
            PriorParams.jitterCovPriorScaleMat(2, 1) = PriorParams.jitterCovPriorScaleMat(1, 2);
        elseif ModelSpec.positionJitterModelIndicator ~= 0,
            error('Invalid positionJitterModelIndicator parameter!')
        end
    end
end