The Genomic Rate of Adaptive Amino-Acid Substitution in Drosophila

Nicolas Bierne and Adam Eyre-Walker

Centre for the Study of Evolution & School of Biological Sciences
University of Sussex
BN19QJ Brighton
UK

Parameters

NP: Number of parameters in the model
(Subscript[P, n]): Number of non-synonymous polymorphisms
(Subscript[P, s]): Number of synonymous polymorphisms
(Subscript[D, n]): Number of non-synonymous substitutions
(Subscript[D, s]): Number of synonymous substitutions
L: Number of codons
Θ: synonymous diversity per codon
λ: synonymous substitution rate per codon
ω: non-synonymous to synonymous diversity ratio
η: number of adaptive non-synonymous substitutions per codon
α: proportion of non-synonymous substitutions which are adaptive (α = ((η L)/((Subscript[D, n]))))

\(TraditionalForm\`\(P\_s\)\&^\) : Poisson distributed with mean	Θ L
\(TraditionalForm\`\(P\_n\)\&^\) : Poisson distributed with mean	ω Θ L
\(TraditionalForm\`\(D\_s\)\&^\) : Poisson distributed with mean	λ L
\(TraditionalForm\`\(D\_n\)\&^\) : Poisson distributed with mean	(ωλ+η)L  or  ((ω λ L)/(1- α))

Data Matrix {{(Subscript[P, n1]),(Subscript[P, s1]),(Subscript[D, n1]),(Subscript[D, s1]),(Subscript[L, 1])},{(Subscript[P, n2]),(Subscript[P, s2]),(Subscript[D, n2]),(Subscript[D, s2]),(Subscript[L, 2])},...,{(Subscript[P, ni]),(Subscript[P, si]),(Subscript[D, ni]),(Subscript[D, si]),(Subscript[L, i])},...}

Dataset

Needs["Statistics`DiscreteDistributions`"]
Needs["Statistics`ContinuousDistributions`"] 
Needs["Statistics`DescriptiveStatistics`"] 
Needs["Statistics`ConfidenceIntervals`"] 
Needs["Graphics`Graphics`"]

D35sim-sim/yak (Smith and Eyre-Walker 2002)

data=(({{0, 11, 9, 28, 256}, {11, 5, 37, 34, 225}, {34, 31, 85, 77, 362}, {17, 12, 165, 80, 227}, {0, 45, 1, 44, 425}, {1, 0, 30, 63, 319}, {7, 6, 34, 57, 406}, {8, 13, 82, 191, 470}, {4, 20, 10, 75, 316}, {12, 22, 116, 112, 534}, {2, 11, 2, 55, 336}, {0, 13, 33, 70, 518}, {4, 2, 44, 78, 335}, {1, 10, 1, 65, 361}, {0, 10, 11, 67, 454}, {1, 26, 43, 144, 458}, {3, 15, 46, 125, 453}, {2, 8, 0, 30, 430}, {1, 10, 12, 38, 215}, {3, 5, 92, 97, 299}, {0, 30, 8, 58, 366}, {8, 22, 104, 104, 454}, {1, 21, 12, 65, 352}, {6, 40, 2, 39, 158}, {2, 15, 5, 47, 302}, {2, 16, 10, 95, 557}, {5, 64, 8, 103, 560}, {8, 50, 21, 100, 453}, {0, 15, 0, 55, 392}, {0, 8, 1, 75, 411}, {1, 12, 0, 31, 246}, {0, 30, 3, 37, 178}, {1, 23, 17, 104, 475}, {0, 2, 2, 3, 348}, {0, 11, 5, 43, 216}}));

Model 1a (NP=2): (Subscript[Θ, i])=Θ, (Subscript[λ, i])=λ, (Subscript[ω, i])=1, (Subscript[α, i])=0

Initialisation

Initialisation of parameters and tables

nbstep = 3000;
TableOfProportion = Array[f, 2];
TableOfTheta = 
        Array[Θ, 2];
TableOfLambda = Array[λ, 2];
TableOfAdaptive = Array[a, 2];
parametro = Array[step, {4, 2}];
Do[
        step[i, j] = 1, {i, 1, 4}, {j, 1, 2}];
proportion = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
theta = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
lambda = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
adaptive = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
metroLogL = Table[0, {j, 1, nbstep}];

Do[proportion[[1, k]] = 0.5; 
        f[k] = proportion[[1, k]];, {k, 1, 
          Length[proportion[[1]]]}];
Do[
        theta[[1, k]] = 0.5; Θ[k] = theta[[1, k]];, {k, 
          1, Length[theta[[1]]]}];
Do[
        lambda[[1, k]] = 0.5; λ[k] = lambda[[1, k]];, {k, 1, 
          Length[lambda[[1]]]}];
Do[adaptive[[1, k]] = 0.5; 
        a[k] = adaptive[[1, k]];, {k, 1, 
          Length[adaptive[[1]]]}];

Loglikelihood function

Loglikelihood := 
Sum[
Log[
PDF[
                PoissonDistribution[Θ[1]*data[[i, 5]]], 
                data[[i, 1]]]*
PDF[
                PoissonDistribution[Θ[1]*data[[i, 5]]], 
                data[[i, 2]]]*
PDF[
                PoissonDistribution[λ[1]*data[[i, 5]]], 
                data[[i, 3]]]*
PDF[
                PoissonDistribution[λ[1]*data[[i, 5]]], 
                data[[i, 4]]]
]
, {i, 
            1, Length[data]}];

metroLogL[[1]] = Loglikelihood;

Metropolis implementation

\(\(implementation[j_, temp_] := 
        Module[
{oldll, newll, oldpara, newpara, threshold, 
            d}, 

threshold[d_] := 
            If[temp <= 0, \((d > 0)\), 
              Exp[d/temp] > 
                Random[]]; 

oldll = 
            metroLogL[\([j - 
                  1]\)]; 

Do[
oldpara = 
              theta[\([j - 1, k]\)]; 
newpara = 
              N[theta[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[2, k]]; 
Θ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((theta[\([j, k]\)] = oldpara; 
                If[step[2, k] > 10\^\(-10\), 
                  step[2, k] *= 0.95]; Θ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[newll - oldll], \((theta[\([j, k]\)] = newpara; 
                      step[2, k] /= 0.95; 
                      oldll = newll; Θ[k] = 
                        newpara)\)
, \((theta[\([j, k]\)] = 
                        oldpara; 
                      If[step[2, k] > 10\^\(-10\), 
                        step[2, k] *= 0.95]; Θ[k] = 
                        oldpara)\)]
, \((theta[\([j, 
                          k]\)] = oldpara; 
                    If[step[2, k] > 10\^\(-10\), 
                      step[2, k] *= 0.95]; Θ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                theta[\([1]\)]]}]; 

Do[
oldpara = 
              lambda[\([j - 1, k]\)]; 
newpara = 
              N[lambda[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[3, k]]; 
λ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((lambda[\([j, k]\)] = oldpara; 
                If[step[3, k] > 10\^\(-10\), 
                  step[3, k] *= 0.95]; λ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((lambda[\([j, k]\)] = newpara; 
                      step[3, k] /= 0.95; 
                      oldll = newll; λ[k] = 
                        newpara)\)
, \((lambda[\([j, 
                            k]\)] = oldpara; 
                      If[step[3, k] > 10\^\(-10\), 
                        step[3, k] *= 0.95]; λ[k] = 
                        oldpara)\)]
, \((lambda[\([j, 
                          k]\)] = oldpara; 
                    If[step[3, k] > 10\^\(-10\), 
                      step[3, k] *= 0.95]; λ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                lambda[\([1]\)]]}]; 

metroLogL[\([j]\)] = 
            oldll; 

If[\((\((j/
                      100)\) == Round[j/100])\) || \((j == 2)\), 
            Print["Log(L) = ", 
              metroLogL[\([j]\)], ";  Θ = ", 
              theta[\([j, 1]\)], ";  λ = ", 
              lambda[\([j, 
                  1]\)]]];

];\)\)

Metropolis

Do[implementation[j, 10], {j, 2, 200}];
Do[implementation[j, 5], {j, 201, 400}];
Do[
        implementation[j, 1], {j, 401, 1400}];
Do[
        implementation[j, 0.8], {j, 1401, 
          1500}];
Do[
        implementation[j, 0.6], {j, 1501, 
          1600}];
Do[
        implementation[j, 0.4], {j, 1601, 
          1700}];
Do[
        implementation[j, 0.2], {j, 1701, 
          1900}];
Do[implementation[j, 0], {j, 1901, nbstep}];

Log(L) = -2037.5;  Θ = 0.03;  λ = 0.14

Model 1b (NP=3): (Subscript[Θ, i])=Θ, (Subscript[λ, i])=λ, (Subscript[ω, i])=ω, (Subscript[α, i])=0

Initialisation

Initialisation of parameters and tables

nbstep = 3000;
TableOfProportion = Array[f, 2];
TableOfTheta = 
        Array[Θ, 2];
TableOfLambda = Array[λ, 2];
TableOfAdaptive = Array[a, 2];
parametro = Array[step, {4, 2}];
Do[
        step[i, j] = 1, {i, 1, 4}, {j, 1, 2}];
proportion = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
theta = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
lambda = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
adaptive = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
metroLogL = Table[0, {j, 1, nbstep}];

Do[proportion[[1, k]] = 0.5; 
        f[k] = proportion[[1, k]];, {k, 1, 
          Length[proportion[[1]]]}];
Do[
        theta[[1, k]] = 0.5; Θ[k] = theta[[1, k]];, {k, 
          1, Length[theta[[1]]]}];
Do[
        lambda[[1, k]] = 0.5; λ[k] = lambda[[1, k]];, {k, 1, 
          Length[lambda[[1]]]}];
Do[adaptive[[1, k]] = 0.5; 
        a[k] = adaptive[[1, k]];, {k, 1, 
          Length[adaptive[[1]]]}];

Loglikelihood function

Loglikelihood := 
Sum[
Log[
PDF[
                PoissonDistribution[
                  f[1]*Θ[1]*data[[i, 5]]], 
                data[[i, 1]]]*
PDF[
                PoissonDistribution[Θ[1]*data[[i, 5]]], 
                data[[i, 2]]]*
PDF[
                PoissonDistribution[f[1]*λ[1]*data[[i, 5]]], 
                data[[i, 3]]]*
PDF[
                PoissonDistribution[λ[1]*data[[i, 5]]], 
                data[[i, 4]]]
]
, {i, 
            1, Length[data]}];

metroLogL[[1]] = Loglikelihood;

Metropolis implementation

\(\(implementation[j_, temp_] := 
        Module[
{oldll, newll, oldpara, newpara, threshold, 
            d}, 

threshold[d_] := 
            If[temp <= 0, \((d > 0)\), 
              Exp[d/temp] > 
                Random[]]; 

oldll = 
            metroLogL[\([j - 
                  1]\)]; 

Do[
oldpara = 
              proportion[\([j - 1, k]\)]; 
newpara = 
              N[proportion[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[1, k]]; 
f[k] = 
              newpara; 
If[
              newpara ≤ 0, \((proportion[\([j, k]\)] = oldpara; 
                If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                f[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((proportion[\([j, k]\)] = newpara; 
                      step[1, k] /= 0.95; oldll = newll; 
                      f[k] = newpara)\)
, \((proportion[\([j, k]\)] = oldpara; 
                      If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                      f[k] = oldpara)\)]
, \((proportion[\([j, k]\)] = oldpara; 
                    If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                    f[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[proportion[\([1]\)]]}]; 

Do[
oldpara = 
              theta[\([j - 1, k]\)]; 
newpara = 
              N[theta[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[2, k]]; 
Θ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((theta[\([j, k]\)] = oldpara; 
                If[step[2, k] > 10\^\(-10\), 
                  step[2, k] *= 0.95]; Θ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[newll - oldll], \((theta[\([j, k]\)] = newpara; 
                      step[2, k] /= 0.95; 
                      oldll = newll; Θ[k] = 
                        newpara)\)
, \((theta[\([j, k]\)] = 
                        oldpara; 
                      If[step[2, k] > 10\^\(-10\), 
                        step[2, k] *= 0.95]; Θ[k] = 
                        oldpara)\)]
, \((theta[\([j, 
                          k]\)] = oldpara; 
                    If[step[2, k] > 10\^\(-10\), 
                      step[2, k] *= 0.95]; Θ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                theta[\([1]\)]]}]; 

Do[
oldpara = 
              lambda[\([j - 1, k]\)]; 
newpara = 
              N[lambda[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[3, k]]; 
λ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((lambda[\([j, k]\)] = oldpara; 
                If[step[3, k] > 10\^\(-10\), 
                  step[3, k] *= 0.95]; λ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((lambda[\([j, k]\)] = newpara; 
                      step[3, k] /= 0.95; 
                      oldll = newll; λ[k] = 
                        newpara)\)
, \((lambda[\([j, 
                            k]\)] = oldpara; 
                      If[step[3, k] > 10\^\(-10\), 
                        step[3, k] *= 0.95]; λ[k] = 
                        oldpara)\)]
, \((lambda[\([j, 
                          k]\)] = oldpara; 
                    If[step[3, k] > 10\^\(-10\), 
                      step[3, k] *= 0.95]; λ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                lambda[\([1]\)]]}]; 

metroLogL[\([j]\)] = 
            oldll; 

If[\((\((j/
                      100)\) == Round[j/100])\) || \((j == 2)\), 
            Print["Log(L) = ", 
              metroLogL[\([j]\)], ";  Θ = ", 
              theta[\([j, 1]\)], ";  λ = ", 
              lambda[\([j, 1]\)], ";  ω = ", 
              proportion[\([j, 
                  1]\)]]];

];\)\)

Metropolis running

Do[implementation[j, 10], {j, 2, 200}];
Do[implementation[j, 5], {j, 201, 400}];
Do[
        implementation[j, 1], {j, 401, 1400}];
Do[
        implementation[j, 0.8], {j, 1401, 
          1500}];
Do[
        implementation[j, 0.6], {j, 1501, 
          1600}];
Do[
        implementation[j, 0.4], {j, 1601, 
          1700}];
Do[
        implementation[j, 0.2], {j, 1701, 
          1900}];
Do[implementation[j, 0], {j, 1901, nbstep}];

Log(L) = -1592.1;  Θ = 0.044;  λ = 0.20;  ω = 0.38

Model 1d (NP=4): (Subscript[Θ, i])=Θ, (Subscript[λ, i])=λ, (Subscript[ω, i])=ω, (Subscript[α, i])=α

Initialisation

Initialisation of parameters and tables

nbstep = 3000;
TableOfProportion = Array[f, 2];
TableOfTheta = 
        Array[Θ, 2];
TableOfLambda = Array[λ, 2];
TableOfAdaptive = Array[a, 2];
parametro = Array[step, {4, 2}];
Do[
        step[i, j] = 1, {i, 1, 4}, {j, 1, 2}];
proportion = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
theta = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
lambda = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
adaptive = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
metroLogL = Table[0, {j, 1, nbstep}];

Do[proportion[[1, k]] = 0.5; 
        f[k] = proportion[[1, k]];, {k, 1, 
          Length[proportion[[1]]]}];
Do[
        theta[[1, k]] = 0.5; Θ[k] = theta[[1, k]];, {k, 
          1, Length[theta[[1]]]}];
Do[
        lambda[[1, k]] = 0.5; λ[k] = lambda[[1, k]];, {k, 1, 
          Length[lambda[[1]]]}];
Do[adaptive[[1, k]] = 0.5; 
        a[k] = adaptive[[1, k]];, {k, 1, 
          Length[adaptive[[1]]]}];

Loglikelihood function

Loglikelihood := 
Sum[
Log[
PDF[
                PoissonDistribution[
                  f[1]*Θ[1]*data[[i, 5]]], 
                data[[i, 1]]]*
PDF[
                PoissonDistribution[Θ[1]*data[[i, 5]]], 
                data[[i, 2]]]*
PDF[
                PoissonDistribution[(f[1]*λ[1]*
                        data[[i, 5]])/(1 - a[1])], 
                data[[i, 3]]]*
PDF[
                PoissonDistribution[λ[1]*data[[i, 5]]], 
                data[[i, 4]]]
]
, {i, 
            1, Length[data]}];

metroLogL[[1]] = Loglikelihood;

Metropolis implementation

\(\(implementation[j_, temp_] := 
        Module[
{oldll, newll, oldpara, newpara, threshold, 
            d}, 

threshold[d_] := 
            If[temp <= 0, \((d > 0)\), 
              Exp[d/temp] > 
                Random[]]; 

oldll = 
            metroLogL[\([j - 
                  1]\)]; 

Do[
oldpara = 
              proportion[\([j - 1, k]\)]; 
newpara = 
              N[proportion[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[1, k]]; 
f[k] = 
              newpara; 
If[
              newpara ≤ 0, \((proportion[\([j, k]\)] = oldpara; 
                If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                f[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((proportion[\([j, k]\)] = newpara; 
                      step[1, k] /= 0.95; oldll = newll; 
                      f[k] = newpara)\)
, \((proportion[\([j, k]\)] = oldpara; 
                      If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                      f[k] = oldpara)\)]
, \((proportion[\([j, k]\)] = oldpara; 
                    If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                    f[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[proportion[\([1]\)]]}]; 

Do[
oldpara = 
              theta[\([j - 1, k]\)]; 
newpara = 
              N[theta[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[2, k]]; 
Θ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((theta[\([j, k]\)] = oldpara; 
                If[step[2, k] > 10\^\(-10\), 
                  step[2, k] *= 0.95]; Θ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[newll - oldll], \((theta[\([j, k]\)] = newpara; 
                      step[2, k] /= 0.95; 
                      oldll = newll; Θ[k] = 
                        newpara)\)
, \((theta[\([j, k]\)] = 
                        oldpara; 
                      If[step[2, k] > 10\^\(-10\), 
                        step[2, k] *= 0.95]; Θ[k] = 
                        oldpara)\)]
, \((theta[\([j, 
                          k]\)] = oldpara; 
                    If[step[2, k] > 10\^\(-10\), 
                      step[2, k] *= 0.95]; Θ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                theta[\([1]\)]]}]; 

Do[
oldpara = 
              lambda[\([j - 1, k]\)]; 
newpara = 
              N[lambda[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[3, k]]; 
λ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((lambda[\([j, k]\)] = oldpara; 
                If[step[3, k] > 10\^\(-10\), 
                  step[3, k] *= 0.95]; λ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((lambda[\([j, k]\)] = newpara; 
                      step[3, k] /= 0.95; 
                      oldll = newll; λ[k] = 
                        newpara)\)
, \((lambda[\([j, 
                            k]\)] = oldpara; 
                      If[step[3, k] > 10\^\(-10\), 
                        step[3, k] *= 0.95]; λ[k] = 
                        oldpara)\)]
, \((lambda[\([j, 
                          k]\)] = oldpara; 
                    If[step[3, k] > 10\^\(-10\), 
                      step[3, k] *= 0.95]; λ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                lambda[\([1]\)]]}]; 

Do[
oldpara = 
              adaptive[\([j - 1, k]\)]; 
newpara = 
              N[adaptive[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[4, k]]; 
a[k] = 
              newpara; 
If[
              newpara ≤ \(-1\), \((adaptive[\([j, k]\)] = oldpara; 
                If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                a[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((adaptive[\([j, k]\)] = newpara; 
                      step[4, k] /= 0.95; oldll = newll; 
                      a[k] = newpara)\)
, \((adaptive[\([j, 
                            k]\)] = oldpara; 
                      If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                      a[k] = oldpara)\)]
, \((adaptive[\([j, k]\)] = oldpara; 
                    If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                    a[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[adaptive[\([1]\)]]}]; 

metroLogL[\([j]\)] = 
            oldll; 

If[\((\((j/
                      100)\) == Round[j/100])\) || \((j == 2)\), 
            Print["Log(L) = ", 
              metroLogL[\([j]\)], ";  Θ = ", 
              theta[\([j, 1]\)], ";  λ = ", 
              lambda[\([j, 1]\)], ";  ω = ", 
              proportion[\([j, 1]\)], ";  α = ", 
              adaptive[\([j, 
                  1]\)]]];

];\)\)

Metropolis running

Do[implementation[j, 10], {j, 2, 200}];
Do[implementation[j, 5], {j, 201, 400}];
Do[
        implementation[j, 1], {j, 401, 1400}];
Do[
        implementation[j, 0.8], {j, 1401, 
          1500}];
Do[
        implementation[j, 0.6], {j, 1501, 
          1600}];
Do[
        implementation[j, 0.4], {j, 1601, 
          1700}];
Do[
        implementation[j, 0.2], {j, 1701, 
          1900}];
Do[implementation[j, 0], {j, 1901, nbstep}];

Log(L) = -1571.2;  Θ = 0.05;  λ = 0.19;  ω = 0.23;  α = 0.46

Model 2a (NP=3n): All (Subscript[Θ, i]), All (Subscript[λ, i]), All (Subscript[ω, i]), α=0

Initialisation

Initialisation of parameters and tables

nbstep = 3000;
TableOfProportion = 
        Array[f, Length[data]];
TableOfTheta = 
        Array[Θ, Length[data]];
TableOfLambda = 
        Array[λ, Length[data]];
parametro = 
        Array[step, {3, Length[data]}];
Do[
        step[i, j] = 1, {i, 1, 3}, {j, 1, 
          Length[data]}];
proportion = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
theta = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
lambda = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
metroLogL = Table[0, {j, 1, nbstep}];

Do[proportion[[1, k]] = 0.5; 
        f[k] = proportion[[1, k]];, {k, 1, 
          Length[proportion[[1]]]}];
Do[
        theta[[1, k]] = 0.5; Θ[k] = theta[[1, k]];, {k, 
          1, Length[theta[[1]]]}];
Do[
        lambda[[1, k]] = 0.5; λ[k] = lambda[[1, k]];, {k, 1, 
          Length[lambda[[1]]]}];

Loglikelihood function

Loglikelihood := 
Sum[
Log[
PDF[
                PoissonDistribution[
                  f[i]*Θ[i]*data[[i, 5]]], 
                data[[i, 1]]]*
PDF[
                PoissonDistribution[Θ[i]*data[[i, 5]]], 
                data[[i, 2]]]*
PDF[
                PoissonDistribution[f[i]*λ[i]*data[[i, 5]]], 
                data[[i, 3]]]*
PDF[
                PoissonDistribution[λ[i]*data[[i, 5]]], 
                data[[i, 4]]]
]
, {i, 
            1, Length[data]}];

metroLogL[[1]] = Loglikelihood;

Metropolis implementation

\(\(implementation[j_, temp_] := 
        Module[
{oldll, newll, oldpara, newpara, threshold, 
            d}, 

threshold[d_] := 
            If[temp <= 0, \((d > 0)\), 
              Exp[d/temp] > 
                Random[]]; 

oldll = 
            metroLogL[\([j - 
                  1]\)]; 

Do[
oldpara = 
              proportion[\([j - 1, k]\)]; 
newpara = 
              N[proportion[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[1, k]]; 
f[k] = 
              newpara; 
If[
              newpara ≤ 0, \((proportion[\([j, k]\)] = oldpara; 
                If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                f[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((proportion[\([j, k]\)] = newpara; 
                      step[1, k] /= 0.95; oldll = newll; 
                      f[k] = newpara)\)
, \((proportion[\([j, k]\)] = oldpara; 
                      If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                      f[k] = oldpara)\)]
, \((proportion[\([j, k]\)] = oldpara; 
                    If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                    f[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[proportion[\([1]\)]]}]; 

Do[
oldpara = 
              theta[\([j - 1, k]\)]; 
newpara = 
              N[theta[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[2, k]]; 
Θ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((theta[\([j, k]\)] = oldpara; 
                If[step[2, k] > 10\^\(-10\), 
                  step[2, k] *= 0.95]; Θ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[newll - oldll], \((theta[\([j, k]\)] = newpara; 
                      step[2, k] /= 0.95; 
                      oldll = newll; Θ[k] = 
                        newpara)\)
, \((theta[\([j, k]\)] = 
                        oldpara; 
                      If[step[2, k] > 10\^\(-10\), 
                        step[2, k] *= 0.95]; Θ[k] = 
                        oldpara)\)]
, \((theta[\([j, 
                          k]\)] = oldpara; 
                    If[step[2, k] > 10\^\(-10\), 
                      step[2, k] *= 0.95]; Θ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                theta[\([1]\)]]}]; 

Do[
oldpara = 
              lambda[\([j - 1, k]\)]; 
newpara = 
              N[lambda[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[3, k]]; 
λ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((lambda[\([j, k]\)] = oldpara; 
                If[step[3, k] > 10\^\(-10\), 
                  step[3, k] *= 0.95]; λ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((lambda[\([j, k]\)] = newpara; 
                      step[3, k] /= 0.95; 
                      oldll = newll; λ[k] = 
                        newpara)\)
, \((lambda[\([j, 
                            k]\)] = oldpara; 
                      If[step[3, k] > 10\^\(-10\), 
                        step[3, k] *= 0.95]; λ[k] = 
                        oldpara)\)]
, \((lambda[\([j, 
                          k]\)] = oldpara; 
                    If[step[3, k] > 10\^\(-10\), 
                      step[3, k] *= 0.95]; λ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                lambda[\([1]\)]]}]; 

metroLogL[\([j]\)] = 
            oldll; 

If[\((\((j/
                      100)\) == Round[j/100])\) || \((j == 2)\), 
            Print["Log(L) = ", 
              metroLogL[\([j]\)]]];

];\)\)

Metropolis running

Do[implementation[j, 10], {j, 2, 200}];
Do[implementation[j, 5], {j, 201, 400}];
Do[
        implementation[j, 1], {j, 401, 1400}];
Do[
        implementation[j, 0.8], {j, 1401, 
          1500}];
Do[
        implementation[j, 0.6], {j, 1501, 
          1600}];
Do[
        implementation[j, 0.4], {j, 1601, 
          1700}];
Do[
        implementation[j, 0.2], {j, 1701, 
          1900}];
Do[implementation[j, 0], {j, 1901, nbstep}];

Log(L) = -331.2

Model 2c (NP=3n+1): All (Subscript[Θ, i]), All (Subscript[λ, i]), All (Subscript[ω, i]), α=cte

Initialisation

Initialisation of parameters and tables

nbstep = 3000;
TableOfProportion = 
        Array[f, Length[data]];
TableOfTheta = 
        Array[Θ, Length[data]];
TableOfLambda = 
        Array[λ, Length[data]];
TableOfAdaptive = Array[a, 2];
parametro = 
        Array[step, {4, Length[data]}];
Do[
        step[i, j] = 1, {i, 1, 4}, {j, 1, 
          Length[data]}];
proportion = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
theta = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
lambda = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
adaptive = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
metroLogL = Table[0, {j, 1, nbstep}];

Do[proportion[[1, k]] = 0.5; 
        f[k] = proportion[[1, k]];, {k, 1, 
          Length[proportion[[1]]]}];
Do[
        theta[[1, k]] = 0.5; Θ[k] = theta[[1, k]];, {k, 
          1, Length[theta[[1]]]}];
Do[
        lambda[[1, k]] = 0.5; λ[k] = lambda[[1, k]];, {k, 1, 
          Length[lambda[[1]]]}];
Do[adaptive[[1, k]] = 0.5; 
        a[k] = adaptive[[1, k]];, {k, 1, 
          Length[adaptive[[1]]]}];

Loglikelihood function

Loglikelihood := 
Sum[
Log[
PDF[
                PoissonDistribution[
                  f[i]*Θ[i]*data[[i, 5]]], 
                data[[i, 1]]]*
PDF[
                PoissonDistribution[Θ[i]*data[[i, 5]]], 
                data[[i, 2]]]*
PDF[
                PoissonDistribution[(f[i]*λ[i]*
                        data[[i, 5]])/(1 - a[1])], 
                data[[i, 3]]]*
PDF[
                PoissonDistribution[λ[i]*data[[i, 5]]], 
                data[[i, 4]]]
]
, {i, 
            1, Length[data]}];

metroLogL[[1]] = Loglikelihood;

Metropolis implementation

\(\(implementation[j_, temp_] := 
        Module[
{oldll, newll, oldpara, newpara, threshold, 
            d}, 

threshold[d_] := 
            If[temp <= 0, \((d > 0)\), 
              Exp[d/temp] > 
                Random[]]; 

oldll = 
            metroLogL[\([j - 
                  1]\)]; 

Do[
oldpara = 
              proportion[\([j - 1, k]\)]; 
newpara = 
              N[proportion[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[1, k]]; 
f[k] = 
              newpara; 
If[
              newpara ≤ 0, \((proportion[\([j, k]\)] = oldpara; 
                If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                f[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((proportion[\([j, k]\)] = newpara; 
                      step[1, k] /= 0.95; oldll = newll; 
                      f[k] = newpara)\)
, \((proportion[\([j, k]\)] = oldpara; 
                      If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                      f[k] = oldpara)\)]
, \((proportion[\([j, k]\)] = oldpara; 
                    If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                    f[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[proportion[\([1]\)]]}]; 

Do[
oldpara = 
              theta[\([j - 1, k]\)]; 
newpara = 
              N[theta[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[2, k]]; 
Θ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((theta[\([j, k]\)] = oldpara; 
                If[step[2, k] > 10\^\(-10\), 
                  step[2, k] *= 0.95]; Θ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[newll - oldll], \((theta[\([j, k]\)] = newpara; 
                      step[2, k] /= 0.95; 
                      oldll = newll; Θ[k] = 
                        newpara)\)
, \((theta[\([j, k]\)] = 
                        oldpara; 
                      If[step[2, k] > 10\^\(-10\), 
                        step[2, k] *= 0.95]; Θ[k] = 
                        oldpara)\)]
, \((theta[\([j, 
                          k]\)] = oldpara; 
                    If[step[2, k] > 10\^\(-10\), 
                      step[2, k] *= 0.95]; Θ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                theta[\([1]\)]]}]; 

Do[
oldpara = 
              lambda[\([j - 1, k]\)]; 
newpara = 
              N[lambda[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[3, k]]; 
λ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((lambda[\([j, k]\)] = oldpara; 
                If[step[3, k] > 10\^\(-10\), 
                  step[3, k] *= 0.95]; λ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((lambda[\([j, k]\)] = newpara; 
                      step[3, k] /= 0.95; 
                      oldll = newll; λ[k] = 
                        newpara)\)
, \((lambda[\([j, 
                            k]\)] = oldpara; 
                      If[step[3, k] > 10\^\(-10\), 
                        step[3, k] *= 0.95]; λ[k] = 
                        oldpara)\)]
, \((lambda[\([j, 
                          k]\)] = oldpara; 
                    If[step[3, k] > 10\^\(-10\), 
                      step[3, k] *= 0.95]; λ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                lambda[\([1]\)]]}]; 

Do[
oldpara = 
              adaptive[\([j - 1, k]\)]; 
newpara = 
              N[adaptive[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[4, k]]; 
a[k] = 
              newpara; 
If[
              newpara ≤ \(-1\), \((adaptive[\([j, k]\)] = oldpara; 
                If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                a[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((adaptive[\([j, k]\)] = newpara; 
                      step[4, k] /= 0.95; oldll = newll; 
                      a[k] = newpara)\)
, \((adaptive[\([j, 
                            k]\)] = oldpara; 
                      If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                      a[k] = oldpara)\)]
, \((adaptive[\([j, k]\)] = oldpara; 
                    If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                    a[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[adaptive[\([1]\)]]}]; 


metroLogL[\([j]\)] = 
            oldll; 

If[\((\((j/
                      100)\) == Round[j/100])\) || \((j == 2)\), 
            Print["Log(L) = ", metroLogL[\([j]\)], ";  α = ", 
              adaptive[\([j, 
                  1]\)]]];

];\)\)

Metropolis running

Do[implementation[j, 10], {j, 2, 200}];
Do[implementation[j, 5], {j, 201, 400}];
Do[
        implementation[j, 1], {j, 401, 1400}];
Do[
        implementation[j, 0.8], {j, 1401, 
          1500}];
Do[
        implementation[j, 0.6], {j, 1501, 
          1600}];
Do[
        implementation[j, 0.4], {j, 1601, 
          1700}];
Do[
        implementation[j, 0.2], {j, 1701, 
          1900}];
Do[implementation[j, 0], {j, 1901, nbstep}];

Log(L) = -327.5;  α = 0.26

Model 2c (NP=3n+1): All (Subscript[Θ, i]), All (Subscript[λ, i]), All (Subscript[ω, i]), α=BetaDistribution(α1,α2)

Initialisation

Initialisation of parameters and tables

nbstep = 2000;
nbclass = 10;
TableOfProportion = 
        Array[f, Length[data]];
TableOfTheta = 
        Array[Θ, Length[data]];
TableOfLambda = 
        Array[λ, Length[data]];
TableOfAdaptive = Array[a, 2];
parametro = 
        Array[step, {4, Length[data]}];
Do[
        step[i, j] = 1, {i, 1, 4}, {j, 1, 
          Length[data]}];
proportion = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
theta = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
lambda = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
adaptive = 
        Table[0, {i, 1, nbstep}, {j, 1, 2}];
metroLogL = Table[0, {j, 1, nbstep}];

Do[proportion[[1, k]] = 0.1; 
        f[k] = proportion[[1, k]];, {k, 1, 
          Length[proportion[[1]]]}];
Do[
        theta[[1, k]] = 0.05; Θ[k] = 
          theta[[1, k]];, {k, 1, 
          Length[theta[[1]]]}];
Do[
        lambda[[1, k]] = 0.2; λ[k] = lambda[[1, k]];, {k, 1, 
          Length[lambda[[1]]]}];
Do[adaptive[[1, k]] = 1; 
        a[k] = adaptive[[1, k]];, {k, 1, 
          Length[adaptive[[1]]]}];

Loglikelihood function

Loglikelihood := 
Sum[
Log[
(1/
                  nbclass)*
Sum[
PDF[
                    PoissonDistribution[
                      f[i]*Θ[i]*data[[i, 5]]], 
                    data[[i, 1]]]*
PDF[
                    PoissonDistribution[Θ[i]*data[[i, 5]]], 
                    data[[i, 2]]]*
PDF[
                    PoissonDistribution[(f[i]*λ[i]*
                            data[[i, 5]])/(1 - 
                            InverseBetaRegularized[((2*aaa + 1)/(2*
                                      nbclass)), a[1], a[2]])], 
                    data[[i, 3]]]*
PDF[
                    PoissonDistribution[λ[i]*data[[i, 5]]], 
                    data[[i, 4]]]
, {aaa, 
                  0, (nbclass - 
                      1)}]
]
, {i, 1, 
            Length[data]}];

metroLogL[[1]] = Loglikelihood;

Metropolis implementation

\(\(implementation[j_, temp_] := 
        Module[
{oldll, newll, oldpara, newpara, threshold, 
            d}, 

threshold[d_] := 
            If[temp <= 0, \((d > 0)\), 
              Exp[d/temp] > 
                Random[]]; 

oldll = 
            metroLogL[\([j - 
                  1]\)]; 

Do[
oldpara = 
              proportion[\([j - 1, k]\)]; 
newpara = 
              N[proportion[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[1, k]]; 
f[k] = 
              newpara; 
If[
              newpara ≤ 0, \((proportion[\([j, k]\)] = oldpara; 
                If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                f[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((proportion[\([j, k]\)] = newpara; 
                      step[1, k] /= 0.95; oldll = newll; 
                      f[k] = newpara)\)
, \((proportion[\([j, k]\)] = oldpara; 
                      If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                      f[k] = oldpara)\)]
, \((proportion[\([j, k]\)] = oldpara; 
                    If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                    f[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[proportion[\([1]\)]]}]; 

Do[
oldpara = 
              theta[\([j - 1, k]\)]; 
newpara = 
              N[theta[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[2, k]]; 
Θ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((theta[\([j, k]\)] = oldpara; 
                If[step[2, k] > 10\^\(-10\), 
                  step[2, k] *= 0.95]; Θ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[newll - oldll], \((theta[\([j, k]\)] = newpara; 
                      step[2, k] /= 0.95; 
                      oldll = newll; Θ[k] = 
                        newpara)\)
, \((theta[\([j, k]\)] = 
                        oldpara; 
                      If[step[2, k] > 10\^\(-10\), 
                        step[2, k] *= 0.95]; Θ[k] = 
                        oldpara)\)]
, \((theta[\([j, 
                          k]\)] = oldpara; 
                    If[step[2, k] > 10\^\(-10\), 
                      step[2, k] *= 0.95]; Θ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                theta[\([1]\)]]}]; 

Do[
oldpara = 
              lambda[\([j - 1, k]\)]; 
newpara = 
              N[lambda[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[3, k]]; 
λ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((lambda[\([j, k]\)] = oldpara; 
                If[step[3, k] > 10\^\(-10\), 
                  step[3, k] *= 0.95]; λ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((lambda[\([j, k]\)] = newpara; 
                      step[3, k] /= 0.95; 
                      oldll = newll; λ[k] = 
                        newpara)\)
, \((lambda[\([j, 
                            k]\)] = oldpara; 
                      If[step[3, k] > 10\^\(-10\), 
                        step[3, k] *= 0.95]; λ[k] = 
                        oldpara)\)]
, \((lambda[\([j, 
                          k]\)] = oldpara; 
                    If[step[3, k] > 10\^\(-10\), 
                      step[3, k] *= 0.95]; λ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                lambda[\([1]\)]]}]; 

Do[
oldpara = 
              adaptive[\([j - 1, k]\)]; 
newpara = 
              N[adaptive[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[4, k]]; 
a[k] = 
              newpara; 
If[
              newpara ≤ 0, \((adaptive[\([j, k]\)] = oldpara; 
                If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                a[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((adaptive[\([j, k]\)] = newpara; 
                      step[4, k] /= 0.95; oldll = newll; 
                      a[k] = newpara)\)
, \((adaptive[\([j, 
                            k]\)] = oldpara; 
                      If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                      a[k] = oldpara)\)]
, \((adaptive[\([j, k]\)] = oldpara; 
                    If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                    a[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[adaptive[\([1]\)]]}]; 


metroLogL[\([j]\)] = 
            oldll; 

If[\((\((j/
                      50)\) == Round[j/50])\) || \((j == 2)\), 
            Print["Log(L) = ", metroLogL[\([j]\)], ";  α1 = ", 
              adaptive[\([j, 1]\)], ";  α2 = ", 
              adaptive[\([j, 
                  2]\)]]];

];\)\)

Metropolis running

Do[implementation[j, 10], {j, 2, 30}];
Do[implementation[j, 5], {j, 31, 60}];
Do[implementation[j, 1], {j, 61, 150}];
Do[
        implementation[j, 0.8], {j, 151, 200}];
Do[
        implementation[j, 0.6], {j, 201, 300}];
Do[
        implementation[j, 0.4], {j, 301, 400}];
Do[
        implementation[j, 0.2], {j, 401, 600}];
Do[implementation[j, 0], {j, 601, nbstep}];

Plot[
      PDF[BetaDistribution[262.925, 762.134], α], {α, 0, 1}, 
      PlotRange -> {{0.1, 0.4}, {0, 35}}]

Model 3 (NP=4n): All (Subscript[Θ, i]), All (Subscript[λ, i]), All (Subscript[ω, i]), All (Subscript[α, i])

Initialisation

Initialisation of parameters and tables

nbstep = 3000;
TableOfProportion = 
        Array[f, Length[data]];
TableOfTheta = 
        Array[Θ, Length[data]];
TableOfLambda = 
        Array[λ, Length[data]];
TableOfAdaptive = 
        Array[a, Length[data]];
parametro = 
        Array[step, {4, Length[data]}];
Do[
        step[i, j] = 1, {i, 1, 4}, {j, 1, 
          Length[data]}];
proportion = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
theta = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
lambda = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
adaptive = 
        Table[0, {i, 1, nbstep}, {j, 1, 
            Length[data]}];
metroLogL = Table[0, {j, 1, nbstep}];

Do[proportion[[1, k]] = 0.5; 
        f[k] = proportion[[1, k]];, {k, 1, 
          Length[proportion[[1]]]}];
Do[
        theta[[1, k]] = 0.5; Θ[k] = theta[[1, k]];, {k, 
          1, Length[theta[[1]]]}];
Do[
        lambda[[1, k]] = 0.5; λ[k] = lambda[[1, k]];, {k, 1, 
          Length[lambda[[1]]]}];
Do[adaptive[[1, k]] = 0.5; 
        a[k] = adaptive[[1, k]];, {k, 1, 
          Length[adaptive[[1]]]}];

Loglikelihood function

Loglikelihood := 
Sum[
Log[
PDF[
                PoissonDistribution[
                  f[i]*Θ[i]*data[[i, 5]]], 
                data[[i, 1]]]*
PDF[
                PoissonDistribution[Θ[i]*data[[i, 5]]], 
                data[[i, 2]]]*
PDF[
                PoissonDistribution[(f[i]*λ[i]*
                        data[[i, 5]])/(1 - a[i])], 
                data[[i, 3]]]*
PDF[
                PoissonDistribution[λ[i]*data[[i, 5]]], 
                data[[i, 4]]]
]
, {i, 
            1, Length[data]}];

metroLogL[[1]] = Loglikelihood;

Metropolis implementation

\(\(implementation[j_, temp_] := 
        Module[
{oldll, newll, oldpara, newpara, threshold, 
            d}, 

threshold[d_] := 
            If[temp <= 0, \((d > 0)\), 
              Exp[d/temp] > 
                Random[]]; 

oldll = 
            metroLogL[\([j - 
                  1]\)]; 

Do[
oldpara = 
              proportion[\([j - 1, k]\)]; 
newpara = 
              N[proportion[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[1, k]]; 
f[k] = 
              newpara; 
If[
              newpara ≤ 0, \((proportion[\([j, k]\)] = oldpara; 
                If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                f[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((proportion[\([j, k]\)] = newpara; 
                      step[1, k] /= 0.95; oldll = newll; 
                      f[k] = newpara)\)
, \((proportion[\([j, k]\)] = oldpara; 
                      If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                      f[k] = oldpara)\)]
, \((proportion[\([j, k]\)] = oldpara; 
                    If[step[1, k] > 10\^\(-10\), step[1, k] *= 0.95]; 
                    f[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[proportion[\([1]\)]]}]; 

Do[
oldpara = 
              theta[\([j - 1, k]\)]; 
newpara = 
              N[theta[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[2, k]]; 
Θ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((theta[\([j, k]\)] = oldpara; 
                If[step[2, k] > 10\^\(-10\), 
                  step[2, k] *= 0.95]; Θ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[newll - oldll], \((theta[\([j, k]\)] = newpara; 
                      step[2, k] /= 0.95; 
                      oldll = newll; Θ[k] = 
                        newpara)\)
, \((theta[\([j, k]\)] = 
                        oldpara; 
                      If[step[2, k] > 10\^\(-10\), 
                        step[2, k] *= 0.95]; Θ[k] = 
                        oldpara)\)]
, \((theta[\([j, 
                          k]\)] = oldpara; 
                    If[step[2, k] > 10\^\(-10\), 
                      step[2, k] *= 0.95]; Θ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                theta[\([1]\)]]}]; 

Do[
oldpara = 
              lambda[\([j - 1, k]\)]; 
newpara = 
              N[lambda[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[3, k]]; 
λ[k] = 
              newpara; 
If[
              newpara ≤ 0, \((lambda[\([j, k]\)] = oldpara; 
                If[step[3, k] > 10\^\(-10\), 
                  step[3, k] *= 0.95]; λ[k] = 
                  oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((lambda[\([j, k]\)] = newpara; 
                      step[3, k] /= 0.95; 
                      oldll = newll; λ[k] = 
                        newpara)\)
, \((lambda[\([j, 
                            k]\)] = oldpara; 
                      If[step[3, k] > 10\^\(-10\), 
                        step[3, k] *= 0.95]; λ[k] = 
                        oldpara)\)]
, \((lambda[\([j, 
                          k]\)] = oldpara; 
                    If[step[3, k] > 10\^\(-10\), 
                      step[3, k] *= 0.95]; λ[k] = 
                      oldpara)\)]
)\)];
, {k, 1, Length[
                lambda[\([1]\)]]}]; 

Do[
oldpara = 
              adaptive[\([j - 1, k]\)]; 
newpara = 
              N[adaptive[\([j - 1, k]\)] + \((2  Random[] - 1)\)*
                    step[4, k]]; 
a[k] = 
              newpara; 
If[
              newpara ≤ 0, \((adaptive[\([j, k]\)] = oldpara; 
                If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                a[k] = oldpara)\)
, \((
newll = Loglikelihood; 
If[
                  Head[Round[newll]] === Integer, 
If[
                    threshold[
                      newll - oldll], \((adaptive[\([j, k]\)] = newpara; 
                      step[4, k] /= 0.95; oldll = newll; 
                      a[k] = newpara)\)
, \((adaptive[\([j, 
                            k]\)] = oldpara; 
                      If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                      a[k] = oldpara)\)]
, \((adaptive[\([j, k]\)] = oldpara; 
                    If[step[4, k] > 10\^\(-10\), step[4, k] *= 0.95]; 
                    a[k] = oldpara)\)]
)\)];
, {k, 1, 
              Length[adaptive[\([1]\)]]}]; 


metroLogL[\([j]\)] = 
            oldll; 

If[\((\((j/
                      100)\) == Round[j/100])\) || \((j == 2)\), 
            Print["Log(L) = ", 
              metroLogL[\([j]\)]]];

];\)\)

Metropolis running

Do[implementation[j, 10], {j, 2, 200}];
Do[implementation[j, 5], {j, 201, 400}];
Do[
        implementation[j, 1], {j, 401, 1400}];
Do[
        implementation[j, 0.8], {j, 1401, 
          1500}];
Do[
        implementation[j, 0.6], {j, 1501, 
          1600}];
Do[
        implementation[j, 0.4], {j, 1601, 
          1700}];
Do[
        implementation[j, 0.2], {j, 1701, 
          1900}];
Do[implementation[j, 0], {j, 1901, nbstep}];

Log(L) = -303