mata mata clear

mata
// Density for minimum of Weibull distributed densities
x = (NULL)
theta1 = (NULL)
theta2 = (NULL)
fweimin = (&tostack(),  &x \
                &enter(),    NULL \
                &enter(),    NULL \
                &fwei(),     &theta1 \
                &swapst(),   NULL \
                &Swei(),     &theta2 \
                &product(),  NULL \
                &rotst(),    NULL \
                &enter(),    NULL \
                &Swei(),     &theta1 \
                &swapst(),   NULL \
                &fwei(),     &theta2 \
                &product(),  NULL \
                &add(),      NULL
                )
x = (0, 2, 8 \
     .7, 1, 6)
theta1 = (-2, 2)
theta2 = (0, .3)
rpnfcn(fweimin)

// Integral of density for minimum of Weibull distributed densities
// MC with antithetic sampling
mata stata set seed 9876

ab = (0, 4 \
      0, 2 \
      2, 4 \
      4, 10)
u = uniform(rows(ab), 30)

intfweimin = (&tostack(), &u \
              &MCatnode(), &ab ) \
              fweimin[2..rows(fweimin), .] \
             (&tostack(), &u \
              &MCatwt(), &ab \
              &intres(), NULL)
rpnfcn(intfweimin)

// Tail integral of density for minimum of Weibull distributed densities
// MC with antithetic sampling
mata stata set seed 4757

a = (0\
     2\
     4)
u = uniform(rows(a), 20)

Sweimin = (&tostack(), &u \
           &tlMCatnode(), &a ) \
          fweimin[2..rows(fweimin), .] \
          (&tostack(), &u \
           &tlMCatwt(), &a \
           &intres(), NULL)
rpnfcn(Sweimin)

// Exact result for above tail integral
Sweimin_exact = (&tostack(), &a \
                 &enter(), NULL \
                 &Swei(), &theta1 \
                 &swapst(), NULL \
                 &Swei(), &theta2 \
                 &product(), NULL)
rpnfcn(Sweimin_exact)

// Forward recurrence density based on minimum of two Weibull distributed
// random variables
a0 = J(1, 1, 0)
u0 = uniform(rows(a0), 20)
frfweimin_exact = (&tostack(), &x ) \
                  Sweimin_exact[2..rows(Sweimin_exact), .] \
                  (&tostack(), &u0 \
                   &tlMCatnode(), &a0 ) \
                  Sweimin_exact[2..rows(Sweimin_exact), .] \
                  (&tostack(), &u0 \
                   &tlMCatwt(), &a0 \
                   &intres(), NULL \
                   &divide(), NULL)
x
rpnfcn(frfweimin_exact)


// Tail integral of forward recurrence density based on minimum
// of two Weibull distributed random variables
// MC with antithetic sampling

a0 = J(3, 1, 0)
u0 = uniform(rows(a0), 10)

frSweimin_exact = (&tostack(), &u \
                   &tlMCatnode(), &a ) \
                  Sweimin_exact[2..rows(Sweimin_exact), .] \
                  (&tostack(), &u \
                   &tlMCatwt(), &a \
                   &intres(), NULL \
                   &tostack(), &u0 \
                   &tlMCatnode(), &a0 ) \
                  Sweimin_exact[2..rows(Sweimin_exact), .] \
                  (&tostack(), &u0 \
                   &tlMCatwt(), &a0 \
                   &intres(), NULL \
                   &divide(), NULL)
a
rpnfcn(frSweimin_exact)
end