function  [X,Y] = func_generateHMM(Qstar, pistar, param, kast, n, lawtype)
% Generate n consecutive HMM observations with specified parameters.

disp('%%%');
disp(['Hidden chain simulation of length ' num2str(n)]);
disp('%%%');

% Generation of the hidden chain.
X=ones(1,n);
cumpi = cumsum(pistar);
u = rand(1,1);
for j=2:kast
    X(:,1) = X(:,1) + (j-1).*(u>=cumpi(j-1)).*(u<cumpi(j));
end
for a=1:n-1
    u=rand(1,1);
    cumQ = cumsum(Qstar(X(:,a),:),2);
    for j=2:kast
         X(:,a+1) = X(:,a+1) + (j-1).*(u>=cumQ(:,j-1)).*(u<cumQ(:,j));
    end
end

% Generation of the observations.
Z = zeros(kast,n);
for x_value = 1:kast
    if strcmp(lawtype(x_value), 'beta')
        % Beta laws.
        param1 = param(x_value,1) * ones(1,n);
        param2 = param(x_value,2) * ones(1,n);
        Z(x_value,:) = betarnd(param1, param2);
    elseif strcmp(lawtype(x_value), 'hist')
        % Piecewise constant densities.
        % Choose which piece will be chosen.
        uniformLaw = rand(1,n);
        labelUnif = rand(1,n);
        piecesNumber = size(param);
        piecesNumber = piecesNumber(2)/3; % Get number of pieces to check.
        % Calculate thresholds (to choose piece with a uniform law)
        threshvalues = zeros(kast, piecesNumber);
        threshvalues(:,1) = (param(:,2)-param(:,1)).*param(:,3);
        for j=2:piecesNumber
            threshvalues(:,j) = threshvalues(:,j-1) + (param(:,3*j-1)-param(:,3*j-2)).*param(:,3*j);
        end
        % Generate observations
        Z(x_value,:) = (labelUnif <= threshvalues(X,1)') .* ( uniformLaw .* (param(x_value,2)-param(x_value,1))' + param(x_value,1)' );
        for piece=2:piecesNumber
            Z(x_value,:) = Z(x_value,:) + (labelUnif > threshvalues(x_value,piece-1)') .* ...
                    (labelUnif <= threshvalues(x_value,piece)') .* ...
                    ( uniformLaw .* (param(x_value,3*piece-1)-param(x_value,3*piece-2))' + param(x_value,3*piece-2)' );
        end
    end
end

Y = Z(kast*(0:length(X)-1) + X);
