\textbf{The lip cancer model}

The number of counties is set to 10.

New values for functional parameters depending on alpha:

\begin{verbatim}
! calculating new value(s) for mu

DO i2 = 1, 10
    mu % one_dim(i2) % new_value = EXP(LOG(expected % one_dim(i2) % &
        value)+alpha % one_dim(1) % new_value*x % one_dim(i2) % &
        value/10+sigma % one_dim(1) % value*b % one_dim(i2) % value)
END DO
\end{verbatim}

Acceptation formula for alpha:

\begin{verbatim}
!  P(alpha') / P(alpha)

p_acc = 1
CALL user_normal_freq(0.0_dp, 100.0_dp, alpha % one_dim(1) % new_value, &
    new_frequency)
CALL user_normal_freq(0.0_dp, 100.0_dp, alpha % one_dim(1) % value, frequency)
p_acc = p_acc * new_frequency / frequency

! P( obs | alpha') / P( obs | alpha)

DO i=1, 10
    ifail=0
    IF (mu % one_dim(i) % new_value > 1.0D+6) THEN
        CALL G01BKF( 1.0D+6, obs % one_dim(i) % value, temp_real, temp_real, &
            new_frequency, ifail)
    ELSE
        CALL G01BKF(mu % one_dim(i) % new_value, obs % one_dim(i) % value, &
            temp_real, temp_real, new_frequency, ifail)
    END IF
    IF (ifail /= 0) THEN
        WRITE (*, *) 'Error when trying to calculate frequency of the poisson &
            &distribution, parameters: mu = ', mu % one_dim(i) % new_value
        STOP
    END IF
    ifail=0
    IF (mu % one_dim(i) % value > 1.0D+6) THEN
        CALL G01BKF( 1.0D+6, obs % one_dim(i) % value, temp_real, temp_real, &
            frequency, ifail)
    ELSE
        CALL G01BKF(mu % one_dim(i) % value, obs % one_dim(i) % value, &
            temp_real, temp_real, frequency, ifail)
    END IF
    IF (ifail /= 0) THEN
        WRITE (*, *) 'Error when trying to calculate frequency of the poisson &
            &distribution, parameters: mu = ', mu % one_dim(i) % value
        STOP
    END IF
    p_acc = p_acc * new_frequency / frequency
END DO

! q(alpha', alpha) / q(alpha, alpha')

transition_real = alpha % one_dim(1) % new_value - alpha % one_dim(1) % value
CALL user_normal_freq(0.0_dp, 10.0_dp, -transition_real, frequency)
CALL user_normal_freq(0.0_dp, 10.0_dp, transition_real, new_frequency)
p_acc = p_acc * frequency / new_frequency
\end{verbatim}

Acceptation code for alpha:

\begin{verbatim}
IF (p_acc > prob(prob_index)) THEN
    alpha % one_dim(1) % value = alpha % one_dim(1) % new_value
    alpha % one_dim(1) % successful_changes = alpha % one_dim(1) % &
        successful_changes + 1
    
    DO i2 = 1, 10
        mu % one_dim(i2) % value = mu % one_dim(i2) % new_value
    END DO
END IF
prob_index = prob_index + 1
IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    CALL G05FAF(0.0_dp, 1.0_dp, SIZE(prob), prob)
    prob_index = 1
END IF
\end{verbatim}

Calculating new values for functional parameters depending on sigma:

\begin{verbatim}
! calculating new value(s) for mu

DO i2 = 1, 10
    mu % one_dim(i2) % new_value = EXP(LOG(expected % one_dim(i2) % &
        value)+alpha % one_dim(1) % value*x % one_dim(i2) % value/10+sigma % &
        one_dim(1) % new_value*b % one_dim(i2) % value)
END DO
\end{verbatim}

Acceptation formula for sigma:

\begin{verbatim}
!  P(sigma') / P(sigma)

p_acc = 1
CALL user_gamma_freq(1.0_dp, 1.0_dp, sigma % one_dim(1) % new_value, &
    new_frequency)
CALL user_gamma_freq(1.0_dp, 1.0_dp, sigma % one_dim(1) % value, frequency)
p_acc = p_acc * new_frequency / frequency

! P( obs | sigma') / P( obs | sigma)

DO i=1, 10
    ifail=0
    IF (mu % one_dim(i) % new_value > 1.0D+6) THEN
        CALL G01BKF( 1.0D+6, obs % one_dim(i) % value, temp_real, temp_real, &
            new_frequency, ifail)
    ELSE
        CALL G01BKF(mu % one_dim(i) % new_value, obs % one_dim(i) % value, &
            temp_real, temp_real, new_frequency, ifail)
    END IF
    IF (ifail /= 0) THEN
        WRITE (*, *) 'Error when trying to calculate frequency of the poisson &
            &distribution, parameters: mu = ', mu % one_dim(i) % new_value
        STOP
    END IF
    ifail=0
    IF (mu % one_dim(i) % value > 1.0D+6) THEN
        CALL G01BKF( 1.0D+6, obs % one_dim(i) % value, temp_real, temp_real, &
            frequency, ifail)
    ELSE
        CALL G01BKF(mu % one_dim(i) % value, obs % one_dim(i) % value, &
            temp_real, temp_real, frequency, ifail)
    END IF
    IF (ifail /= 0) THEN
        WRITE (*, *) 'Error when trying to calculate frequency of the poisson &
            &distribution, parameters: mu = ', mu % one_dim(i) % value
        STOP
    END IF
    p_acc = p_acc * new_frequency / frequency
END DO

! q(sigma', sigma) / q(sigma, sigma')

IF (sigma % one_dim(1) % value < 0.0010_dp .OR. sigma % one_dim(1) % value > &
    5.0_dp) THEN
    frequency = 0
ELSE
    frequency = 1.0 / (5.0_dp - (0.0010_dp))
END IF
IF (sigma % one_dim(1) % new_value < 0.0010_dp .OR. sigma % one_dim(1) % &
    new_value > 5.0_dp) THEN
    new_frequency = 0
ELSE
    new_frequency = 1.0 / (5.0_dp - (0.0010_dp))
END IF
p_acc = p_acc * frequency / new_frequency
\end{verbatim}

Acceptation code for sigma (if the update strategy is sequential):

\begin{verbatim}
IF (p_acc > prob(prob_index)) THEN
    sigma % one_dim(1) % value = sigma % one_dim(1) % new_value
    sigma % one_dim(1) % successful_changes = sigma % one_dim(1) % &
        successful_changes + 1
    
    DO i2 = 1, 10
        mu % one_dim(i2) % value = mu % one_dim(i2) % new_value
    END DO
END IF
prob_index = prob_index + 1
IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    CALL G05FAF(0.0_dp, 1.0_dp, SIZE(prob), prob)
    prob_index = 1
END IF
\end{verbatim}

Calculating new values for functional parameters depending on b:

Note that the q:s and t:s of neighbour squares don't need to be updated at this phase. They are updated only if the new value for b is accepted.

\begin{verbatim}
! calculating new value(s) for mu

i2 = i
mu % one_dim(i2) % new_value = EXP(LOG(expected % one_dim(i2) % value)+alpha &
    % one_dim(1) % value*x % one_dim(i2) % value/10+sigma % one_dim(1) % &
    value*b % one_dim(i2) % new_value)
\end{verbatim}

Acceptation formula for b:

\begin{verbatim}
!  P(b') / P(b)

p_acc = 1
CALL user_normal_freq(t % one_dim(i) % value, n % one_dim(i) % value, b % &
    one_dim(i) % new_value, new_frequency)
CALL user_normal_freq(t % one_dim(i) % value, n % one_dim(i) % value, b % &
    one_dim(i) % value, frequency)
p_acc = p_acc * new_frequency / frequency


! P( obs | b') / P( obs | b)

ifail=0
IF (mu % one_dim(i) % new_value > 1.0D+6) THEN
    CALL G01BKF( 1.0D+6, obs % one_dim(i) % value, temp_real, temp_real, &
        new_frequency, ifail)
ELSE
    CALL G01BKF(mu % one_dim(i) % new_value, obs % one_dim(i) % value, &
        temp_real, temp_real, new_frequency, ifail)
END IF
IF (ifail /= 0) THEN
    WRITE (*, *) 'Error when trying to calculate frequency of the poisson &
        &distribution, parameters: mu = ', mu % one_dim(i) % new_value
    STOP
END IF
ifail=0
IF (mu % one_dim(i) % value > 1.0D+6) THEN
    CALL G01BKF( 1.0D+6, obs % one_dim(i) % value, temp_real, temp_real, &
        frequency, ifail)
ELSE
    CALL G01BKF(mu % one_dim(i) % value, obs % one_dim(i) % value, temp_real, &
        temp_real, frequency, ifail)
END IF
IF (ifail /= 0) THEN
    WRITE (*, *) 'Error when trying to calculate frequency of the poisson &
        &distribution, parameters: mu = ', mu % one_dim(i) % value
    STOP
END IF
p_acc = p_acc * new_frequency / frequency

! q(b', b) / q(b, b')

IF (b % one_dim(i) % value < -10.0_dp .OR. b % one_dim(i) % value > 10.0_dp) &
    THEN
    frequency = 0
ELSE
    frequency = 1.0 / (10.0_dp - (-10.0_dp))
END IF
IF (b % one_dim(i) % new_value < -10.0_dp .OR. b % one_dim(i) % new_value > &
    10.0_dp) THEN
    new_frequency = 0
ELSE
    new_frequency = 1.0 / (10.0_dp - (-10.0_dp))
END IF
p_acc = p_acc * frequency / new_frequency
\end{verbatim}

Acceptation code for b (if the update strategy is sequential):

\begin{verbatim}
IF (p_acc > prob(prob_index)) THEN
    b % one_dim(i) % value = b % one_dim(i) % new_value
    b % one_dim(i) % successful_changes = b % one_dim(i) % successful_changes &
        + 1
    
    i2 = i
    mu % one_dim(i2) % value = mu % one_dim(i2) % new_value
    ! 1 to number of neighbours
    DO k = 1, spatial(i, 1)
        sum_real = 0
        neighbour = spatial(i, k + 1)
        ! updating q % one_dim(spatial(i, k + 1))
        DO h = 1, spatial(neighbour, 1)
            sum_real = sum_real + b % one_dim(spatial(neighbour, h + 1)) % &
                value
        END DO
        q % one_dim(neighbour) % value = sum_real
    END DO
    DO k = 1, spatial(i, 1)
        neighbour = spatial(i,  k +1)
        t % one_dim(neighbour) % value = q % one_dim(neighbour) % value/n % &
            one_dim(neighbour) % value
    END DO
END IF
prob_index = prob_index + 1
IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    CALL G05FAF(0.0_dp, 1.0_dp, SIZE(prob), prob)
    prob_index = 1
END IF
\end{verbatim}