diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,4 +1,17 @@
 
+Version 4.2
+-----
+
+* Added new random distributions: lognormal, Gamma, Beta, Weibull and 
+  a generic discrete by pdf.
+
+* The items can be removed from the queue; moreover, the queue can be 
+  cleared.
+
+* Added a simplified API for accessing the results of simulation.
+
+* Added the Gate entity.
+
 Version 4.1.1
 -----
 
diff --git a/Simulation/Aivika.hs b/Simulation/Aivika.hs
--- a/Simulation/Aivika.hs
+++ b/Simulation/Aivika.hs
@@ -23,6 +23,7 @@
         module Simulation.Aivika.Dynamics.Memo.Unboxed,
         module Simulation.Aivika.Dynamics.Random,
         module Simulation.Aivika.Event,
+        module Simulation.Aivika.Gate,
         module Simulation.Aivika.Generator,
         module Simulation.Aivika.Net,
         module Simulation.Aivika.Net.Random,
@@ -65,6 +66,7 @@
 import Simulation.Aivika.Dynamics.Memo.Unboxed
 import Simulation.Aivika.Dynamics.Random
 import Simulation.Aivika.Event
+import Simulation.Aivika.Gate
 import Simulation.Aivika.Generator
 import Simulation.Aivika.Net
 import Simulation.Aivika.Net.Random
diff --git a/Simulation/Aivika/Activity/Random.hs b/Simulation/Aivika/Activity/Random.hs
--- a/Simulation/Aivika/Activity/Random.hs
+++ b/Simulation/Aivika/Activity/Random.hs
@@ -15,19 +15,32 @@
 module Simulation.Aivika.Activity.Random
        (newRandomUniformActivity,
         newRandomUniformIntActivity,
+        newRandomTriangularActivity,
         newRandomNormalActivity,
+        newRandomLogNormalActivity,
         newRandomExponentialActivity,
         newRandomErlangActivity,
         newRandomPoissonActivity,
         newRandomBinomialActivity,
+        newRandomGammaActivity,
+        newRandomBetaActivity,
+        newRandomWeibullActivity,
+        newRandomDiscreteActivity,
         newPreemptibleRandomUniformActivity,
         newPreemptibleRandomUniformIntActivity,
+        newPreemptibleRandomTriangularActivity,
         newPreemptibleRandomNormalActivity,
+        newPreemptibleRandomLogNormalActivity,
         newPreemptibleRandomExponentialActivity,
         newPreemptibleRandomErlangActivity,
         newPreemptibleRandomPoissonActivity,
-        newPreemptibleRandomBinomialActivity) where
+        newPreemptibleRandomBinomialActivity,
+        newPreemptibleRandomGammaActivity,
+        newPreemptibleRandomBetaActivity,
+        newPreemptibleRandomWeibullActivity,
+        newPreemptibleRandomDiscreteActivity) where
 
+import Simulation.Aivika.Generator
 import Simulation.Aivika.Simulation
 import Simulation.Aivika.Process
 import Simulation.Aivika.Process.Random
@@ -62,6 +75,22 @@
   newPreemptibleRandomUniformIntActivity False
 
 -- | Create a new activity that holds the process for a random time interval
+-- having the triangular distribution, when processing every input element.
+--
+-- By default, it is assumed that the activity process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomTriangularActivity :: Double
+                               -- ^ the minimum time interval
+                               -> Double
+                               -- ^ the median of the time interval
+                               -> Double
+                               -- ^ the maximum time interval
+                               -> Simulation (Activity () a a)
+newRandomTriangularActivity =
+  newPreemptibleRandomTriangularActivity False
+
+-- | Create a new activity that holds the process for a random time interval
 -- distributed normally, when processing every input element.
 --
 -- By default, it is assumed that the activity process cannot be preempted,
@@ -76,6 +105,22 @@
   newPreemptibleRandomNormalActivity False
          
 -- | Create a new activity that holds the process for a random time interval
+-- having the lognormal distribution, when processing every input element.
+--
+-- By default, it is assumed that the activity process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomLogNormalActivity :: Double
+                              -- ^ the mean of a normal distribution which
+                              -- this distribution is derived from
+                              -> Double
+                              -- ^ the deviation of a normal distribution which
+                              -- this distribution is derived from
+                              -> Simulation (Activity () a a)
+newRandomLogNormalActivity =
+  newPreemptibleRandomLogNormalActivity False
+
+-- | Create a new activity that holds the process for a random time interval
 -- distributed exponentially with the specified mean (the reciprocal of the rate),
 -- when processing every input element.
 --
@@ -132,6 +177,63 @@
   newPreemptibleRandomBinomialActivity False
 
 -- | Create a new activity that holds the process for a random time interval
+-- having the Gamma distribution with the specified shape and scale,
+-- when processing every input element.
+--
+-- By default, it is assumed that the activity process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomGammaActivity :: Double
+                          -- ^ the shape
+                          -> Double
+                          -- ^ the scale (a reciprocal of the rate)
+                          -> Simulation (Activity () a a)
+newRandomGammaActivity =
+  newPreemptibleRandomGammaActivity False
+
+-- | Create a new activity that holds the process for a random time interval
+-- having the Beta distribution with the specified shape parameters (alpha and beta),
+-- when processing every input element.
+--
+-- By default, it is assumed that the activity process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomBetaActivity :: Double
+                         -- ^ shape (alpha)
+                         -> Double
+                         -- ^ shape (beta)
+                         -> Simulation (Activity () a a)
+newRandomBetaActivity =
+  newPreemptibleRandomBetaActivity False
+
+-- | Create a new activity that holds the process for a random time interval
+-- having the Weibull distribution with the specified shape and scale,
+-- when processing every input element.
+--
+-- By default, it is assumed that the activity process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomWeibullActivity :: Double
+                            -- ^ shape
+                            -> Double
+                            -- ^ scale
+                            -> Simulation (Activity () a a)
+newRandomWeibullActivity =
+  newPreemptibleRandomWeibullActivity False
+
+-- | Create a new activity that holds the process for a random time interval
+-- having the specified discrete distribution, when processing every input element.
+--
+-- By default, it is assumed that the activity process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomDiscreteActivity :: DiscretePDF Double
+                             -- ^ the discrete probability density function
+                             -> Simulation (Activity () a a)
+newRandomDiscreteActivity =
+  newPreemptibleRandomDiscreteActivity False
+
+-- | Create a new activity that holds the process for a random time interval
 -- distributed uniformly, when processing every input element.
 newPreemptibleRandomUniformActivity :: Bool
                                        -- ^ whether the activity process can be preempted
@@ -160,6 +262,22 @@
      return a
 
 -- | Create a new activity that holds the process for a random time interval
+-- having the triangular distribution, when processing every input element.
+newPreemptibleRandomTriangularActivity :: Bool
+                                          -- ^ whether the activity process can be preempted
+                                          -> Double
+                                          -- ^ the minimum time interval
+                                          -> Double
+                                          -- ^ the median of the time interval
+                                          -> Double
+                                          -- ^ the maximum time interval
+                                          -> Simulation (Activity () a a)
+newPreemptibleRandomTriangularActivity preemptible min median max =
+  newPreemptibleActivity preemptible $ \a ->
+  do randomTriangularProcess_ min median max
+     return a
+
+-- | Create a new activity that holds the process for a random time interval
 -- distributed normally, when processing every input element.
 newPreemptibleRandomNormalActivity :: Bool
                                       -- ^ whether the activity process can be preempted
@@ -172,8 +290,24 @@
   newPreemptibleActivity preemptible $ \a ->
   do randomNormalProcess_ mu nu
      return a
-         
+
 -- | Create a new activity that holds the process for a random time interval
+-- having the lognormal distribution, when processing every input element.
+newPreemptibleRandomLogNormalActivity :: Bool
+                                         -- ^ whether the activity process can be preempted
+                                         -> Double
+                                         -- ^ the mean of a normal distribution which
+                                         -- this distribution is derived from
+                                         -> Double
+                                         -- ^ the deviation of a normal distribution which
+                                         -- this distribution is derived from
+                                         -> Simulation (Activity () a a)
+newPreemptibleRandomLogNormalActivity preemptible mu nu =
+  newPreemptibleActivity preemptible $ \a ->
+  do randomLogNormalProcess_ mu nu
+     return a
+
+-- | Create a new activity that holds the process for a random time interval
 -- distributed exponentially with the specified mean (the reciprocal of the rate),
 -- when processing every input element.
 newPreemptibleRandomExponentialActivity :: Bool
@@ -227,4 +361,61 @@
 newPreemptibleRandomBinomialActivity preemptible prob trials =
   newPreemptibleActivity preemptible $ \a ->
   do randomBinomialProcess_ prob trials
+     return a
+
+-- | Create a new activity that holds the process for a random time interval
+-- having the Gamma distribution with the specified shape and scale,
+-- when processing every input element.
+newPreemptibleRandomGammaActivity :: Bool
+                                     -- ^ whether the activity process can be preempted
+                                     -> Double
+                                     -- ^ the shape
+                                     -> Double
+                                     -- ^ the scale
+                                     -> Simulation (Activity () a a)
+newPreemptibleRandomGammaActivity preemptible kappa theta =
+  newPreemptibleActivity preemptible $ \a ->
+  do randomGammaProcess_ kappa theta
+     return a
+
+-- | Create a new activity that holds the process for a random time interval
+-- having the Beta distribution with the specified shape parameters (alpha and beta),
+-- when processing every input element.
+newPreemptibleRandomBetaActivity :: Bool
+                                    -- ^ whether the activity process can be preempted
+                                    -> Double
+                                    -- ^ shape (alpha)
+                                    -> Double
+                                    -- ^ shape (beta)
+                                    -> Simulation (Activity () a a)
+newPreemptibleRandomBetaActivity preemptible alpha beta =
+  newPreemptibleActivity preemptible $ \a ->
+  do randomBetaProcess_ alpha beta
+     return a
+
+-- | Create a new activity that holds the process for a random time interval
+-- having the Weibull distribution with the specified shape and scale,
+-- when processing every input element.
+newPreemptibleRandomWeibullActivity :: Bool
+                                       -- ^ whether the activity process can be preempted
+                                       -> Double
+                                       -- ^ shape
+                                       -> Double
+                                       -- ^ scale
+                                       -> Simulation (Activity () a a)
+newPreemptibleRandomWeibullActivity preemptible alpha beta =
+  newPreemptibleActivity preemptible $ \a ->
+  do randomWeibullProcess_ alpha beta
+     return a
+
+-- | Create a new activity that holds the process for a random time interval
+-- having the specified discrete distribution, when processing every input element.
+newPreemptibleRandomDiscreteActivity :: Bool
+                                        -- ^ whether the activity process can be preempted
+                                        -> DiscretePDF Double
+                                        -- ^ the discrete probability density function
+                                        -> Simulation (Activity () a a)
+newPreemptibleRandomDiscreteActivity preemptible dpdf =
+  newPreemptibleActivity preemptible $ \a ->
+  do randomDiscreteProcess_ dpdf
      return a
diff --git a/Simulation/Aivika/DoubleLinkedList.hs b/Simulation/Aivika/DoubleLinkedList.hs
--- a/Simulation/Aivika/DoubleLinkedList.hs
+++ b/Simulation/Aivika/DoubleLinkedList.hs
@@ -18,10 +18,14 @@
         listAddLast,
         listRemoveFirst,
         listRemoveLast,
+        listRemove,
+        listRemoveBy,
         listFirst,
         listLast) where 
 
 import Data.IORef
+import Data.Maybe
+
 import Control.Monad
 
 -- | A cell of the double-linked list.
@@ -62,7 +66,8 @@
 listInsertFirst :: DoubleLinkedList a -> a -> IO ()
 listInsertFirst x v =
   do size <- readIORef (listSize x)
-     writeIORef (listSize x) (size + 1)
+     let size' = size + 1
+     size' `seq` writeIORef (listSize x) size'
      head <- readIORef (listHead x)
      case head of
        Nothing ->
@@ -86,7 +91,8 @@
 listAddLast :: DoubleLinkedList a -> a -> IO ()
 listAddLast x v =
   do size <- readIORef (listSize x)
-     writeIORef (listSize x) (size + 1)
+     let size' = size + 1
+     size' `seq` writeIORef (listSize x) size'
      tail <- readIORef (listTail x)
      case tail of
        Nothing ->
@@ -115,7 +121,8 @@
          error "Empty list: listRemoveFirst"
        Just h ->
          do size  <- readIORef (listSize x)
-            writeIORef (listSize x) (size - 1)
+            let size' = size - 1
+            size' `seq` writeIORef (listSize x) size'
             head' <- readIORef (itemNext h)
             case head' of
               Nothing ->
@@ -134,7 +141,8 @@
          error "Empty list: listRemoveLast"
        Just t ->
          do size  <- readIORef (listSize x)
-            writeIORef (listSize x) (size - 1)
+            let size' =  size - 1
+            size' `seq` writeIORef (listSize x) size'
             tail' <- readIORef (itemPrev t)
             case tail' of
               Nothing ->
@@ -163,3 +171,39 @@
          error "Empty list: listLast"
        Just t ->
          return $ itemVal t
+
+-- | Remove the specified element from the list and return a flag
+-- indicating whether the element was found and removed.
+listRemove :: Eq a => DoubleLinkedList a -> a -> IO Bool
+listRemove x v = fmap isJust $ listRemoveBy x (== v)
+
+-- | Remove an element satisfying the specified predicate and return
+-- the element if found.
+listRemoveBy :: DoubleLinkedList a -> (a -> Bool) -> IO (Maybe a)
+listRemoveBy x p = readIORef (listHead x) >>= loop
+  where loop item =
+          case item of
+            Nothing   -> return Nothing
+            Just item ->
+              do let f = p (itemVal item)
+                 if not f
+                   then readIORef (itemNext item) >>= loop
+                   else do size <- readIORef (listSize x)
+                           prev <- readIORef (itemPrev item)
+                           next <- readIORef (itemNext item)
+                           let size' = size - 1
+                           size' `seq` writeIORef (listSize x) size'
+                           case (prev, next) of
+                             (Nothing, Nothing) ->
+                               do writeIORef (listHead x) Nothing
+                                  writeIORef (listTail x) Nothing
+                             (Nothing, head' @ (Just item')) ->
+                               do writeIORef (itemPrev item') Nothing
+                                  writeIORef (listHead x) head'
+                             (tail' @ (Just item'), Nothing) ->
+                               do writeIORef (itemNext item') Nothing
+                                  writeIORef (listTail x) tail'
+                             (Just prev', Just next') ->
+                               do writeIORef (itemNext prev') (Just next')
+                                  writeIORef (itemPrev next') (Just prev')
+                           return (Just $ itemVal item)
diff --git a/Simulation/Aivika/Dynamics/Random.hs b/Simulation/Aivika/Dynamics/Random.hs
--- a/Simulation/Aivika/Dynamics/Random.hs
+++ b/Simulation/Aivika/Dynamics/Random.hs
@@ -22,11 +22,17 @@
 module Simulation.Aivika.Dynamics.Random
        (memoRandomUniformDynamics,
         memoRandomUniformIntDynamics,
+        memoRandomTriangularDynamics,
         memoRandomNormalDynamics,
+        memoRandomLogNormalDynamics,
         memoRandomExponentialDynamics,
         memoRandomErlangDynamics,
         memoRandomPoissonDynamics,
-        memoRandomBinomialDynamics) where
+        memoRandomBinomialDynamics,
+        memoRandomGammaDynamics,
+        memoRandomBetaDynamics,
+        memoRandomWeibullDynamics,
+        memoRandomDiscreteDynamics) where
 
 import System.Random
 
@@ -38,9 +44,10 @@
 import Simulation.Aivika.Internal.Simulation
 import Simulation.Aivika.Internal.Dynamics
 import Simulation.Aivika.Dynamics.Memo.Unboxed
+import Simulation.Aivika.Unboxed
 
 -- | Computation that generates random numbers distributed uniformly and
--- memoizes them in the integration time points.
+-- memoizes the numbers in the integration time points.
 memoRandomUniformDynamics :: Dynamics Double     -- ^ minimum
                              -> Dynamics Double  -- ^ maximum
                              -> Simulation (Dynamics Double)
@@ -53,7 +60,7 @@
      generateUniform g min' max'
 
 -- | Computation that generates random integer numbers distributed uniformly and
--- memoizes them in the integration time points.
+-- memoizes the numbers in the integration time points.
 memoRandomUniformIntDynamics :: Dynamics Int     -- ^ minimum
                                 -> Dynamics Int  -- ^ maximum
                                 -> Simulation (Dynamics Int)
@@ -65,8 +72,23 @@
      max' <- invokeDynamics p max
      generateUniformInt g min' max'
 
+-- | Computation that generates random numbers from the triangular distribution
+-- and memoizes the numbers in the integration time points.
+memoRandomTriangularDynamics :: Dynamics Double     -- ^ minimum
+                                -> Dynamics Double  -- ^ median
+                                -> Dynamics Double  -- ^ maximum
+                                -> Simulation (Dynamics Double)
+memoRandomTriangularDynamics min median max =
+  memo0Dynamics $
+  Dynamics $ \p ->
+  do let g = runGenerator $ pointRun p
+     min' <- invokeDynamics p min
+     median' <- invokeDynamics p median
+     max' <- invokeDynamics p max
+     generateTriangular g min' median' max'
+
 -- | Computation that generates random numbers distributed normally and
--- memoizes them in the integration time points.
+-- memoizes the numbers in the integration time points.
 memoRandomNormalDynamics :: Dynamics Double     -- ^ mean
                             -> Dynamics Double  -- ^ deviation
                             -> Simulation (Dynamics Double)
@@ -78,10 +100,27 @@
      nu' <- invokeDynamics p nu
      generateNormal g mu' nu'
 
+-- | Computation that generates random numbers from the lognormal distribution
+-- and memoizes the numbers in the integration time points.
+memoRandomLogNormalDynamics :: Dynamics Double
+                               -- ^ the mean of a normal distribution which
+                               -- this distribution is derived from
+                               -> Dynamics Double
+                               -- ^ the deviation of a normal distribution which
+                               -- this distribution is derived from
+                               -> Simulation (Dynamics Double)
+memoRandomLogNormalDynamics mu nu =
+  memo0Dynamics $
+  Dynamics $ \p ->
+  do let g = runGenerator $ pointRun p
+     mu' <- invokeDynamics p mu
+     nu' <- invokeDynamics p nu
+     generateLogNormal g mu' nu'
+
 -- | Computation that generates exponential random numbers with the specified mean
--- (the reciprocal of the rate) and memoizes them in the integration time points.
+-- (the reciprocal of the rate) and memoizes the numbers in the integration time points.
 memoRandomExponentialDynamics :: Dynamics Double
-                                 -- ^ the mean (the reciprocal of the rate)
+                                 -- ^ the mean (a reciprocal of the rate)
                                  -> Simulation (Dynamics Double)
 memoRandomExponentialDynamics mu =
   memo0Dynamics $
@@ -91,10 +130,10 @@
      generateExponential g mu'
 
 -- | Computation that generates the Erlang random numbers with the specified scale
--- (the reciprocal of the rate) and integer shape but memoizes them in the integration
--- time points.
+-- (the reciprocal of the rate) and integer shape but memoizes the numbers in
+-- the integration time points.
 memoRandomErlangDynamics :: Dynamics Double
-                            -- ^ the scale (the reciprocal of the rate)
+                            -- ^ the scale (a reciprocal of the rate)
                             -> Dynamics Int
                             -- ^ the shape
                             -> Simulation (Dynamics Double)
@@ -107,7 +146,7 @@
      generateErlang g beta' m'
 
 -- | Computation that generats the Poisson random numbers with the specified mean
--- and memoizes them in the integration time points.
+-- and memoizes the numbers in the integration time points.
 memoRandomPoissonDynamics :: Dynamics Double
                              -- ^ the mean
                              -> Simulation (Dynamics Int)
@@ -119,7 +158,7 @@
      generatePoisson g mu'
 
 -- | Computation that generates binomial random numbers with the specified
--- probability and trials but memoizes them in the integration time points.
+-- probability and trials but memoizes the numbers in the integration time points.
 memoRandomBinomialDynamics :: Dynamics Double  -- ^ the probability
                               -> Dynamics Int  -- ^ the number of trials
                               -> Simulation (Dynamics Int)
@@ -130,3 +169,55 @@
      prob' <- invokeDynamics p prob
      trials' <- invokeDynamics p trials
      generateBinomial g prob' trials'
+
+-- | Computation that generates random numbers from the Gamma distribution
+-- with the specified shape and scale but memoizes the numbers in
+-- the integration time points.
+memoRandomGammaDynamics :: Dynamics Double     -- ^ shape
+                           -> Dynamics Double  -- ^ scale (a reciprocal of the rate)
+                           -> Simulation (Dynamics Double)
+memoRandomGammaDynamics kappa theta =
+  memo0Dynamics $
+  Dynamics $ \p ->
+  do let g = runGenerator $ pointRun p
+     kappa' <- invokeDynamics p kappa
+     theta' <- invokeDynamics p theta
+     generateGamma g kappa' theta'
+
+-- | Computation that generates random numbers from the Beta distribution
+-- by the specified shape parameters and memoizes the numbers in
+-- the integration time points.
+memoRandomBetaDynamics :: Dynamics Double     -- ^ shape (alpha)
+                          -> Dynamics Double  -- ^ shape (beta)
+                          -> Simulation (Dynamics Double)
+memoRandomBetaDynamics alpha beta =
+  memo0Dynamics $
+  Dynamics $ \p ->
+  do let g = runGenerator $ pointRun p
+     alpha' <- invokeDynamics p alpha
+     beta'  <- invokeDynamics p beta
+     generateBeta g alpha' beta'
+
+-- | Computation that generates random numbers from the Weibull distribution
+-- with the specified shape and scale but memoizes the numbers in
+-- the integration time points.
+memoRandomWeibullDynamics :: Dynamics Double     -- ^ shape
+                             -> Dynamics Double  -- ^ scale
+                             -> Simulation (Dynamics Double)
+memoRandomWeibullDynamics alpha beta =
+  memo0Dynamics $
+  Dynamics $ \p ->
+  do let g = runGenerator $ pointRun p
+     alpha' <- invokeDynamics p alpha
+     beta'  <- invokeDynamics p beta
+     generateWeibull g alpha' beta'
+
+-- | Computation that generates random values from the specified discrete
+-- distribution and memoizes the values in the integration time points.
+memoRandomDiscreteDynamics :: Unboxed a => Dynamics (DiscretePDF a) -> Simulation (Dynamics a)
+memoRandomDiscreteDynamics dpdf =
+  memo0Dynamics $
+  Dynamics $ \p ->
+  do let g = runGenerator $ pointRun p
+     dpdf' <- invokeDynamics p dpdf
+     generateDiscrete g dpdf'
diff --git a/Simulation/Aivika/Gate.hs b/Simulation/Aivika/Gate.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Gate.hs
@@ -0,0 +1,91 @@
+
+-- |
+-- Module     : Simulation.Aivika.Gate
+-- Copyright  : Copyright (c) 2009-2015, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.1
+--
+-- The module defines a gate which can be either opened or closed.
+--
+module Simulation.Aivika.Gate
+       (Gate,
+        newGate,
+        newGateOpened,
+        newGateClosed,
+        openGate,
+        closeGate,
+        gateOpened,
+        gateClosed,
+        awaitGateOpened,
+        awaitGateClosed,
+        gateChanged_) where
+
+import Control.Monad
+
+import Simulation.Aivika.Simulation
+import Simulation.Aivika.Event
+import Simulation.Aivika.Process
+import Simulation.Aivika.Signal
+import Simulation.Aivika.Ref
+
+-- | Represents a gate, which can be either opened or closed.
+data Gate = Gate { gateRef :: Ref Bool }
+
+-- | Create a new gate, specifying whether the gate is initially open.
+newGate :: Bool -> Simulation Gate
+newGate opened =
+  do r <- newRef opened
+     return Gate { gateRef = r }
+
+-- | Create a new initially open gate.
+newGateOpened :: Simulation Gate
+newGateOpened = newGate True
+
+-- | Create a new initially close gate.
+newGateClosed :: Simulation Gate
+newGateClosed = newGate False
+
+-- | Open the gate if it was closed.
+openGate :: Gate -> Event ()
+openGate gate =
+  writeRef (gateRef gate) True
+
+-- | Close the gate if it was open.
+closeGate :: Gate -> Event ()
+closeGate gate =
+  writeRef (gateRef gate) False
+
+-- | Test whether the gate is open.
+gateOpened :: Gate -> Event Bool
+gateOpened gate =
+  readRef (gateRef gate)
+
+-- | Test whether the gate is closed.
+gateClosed :: Gate -> Event Bool
+gateClosed gate =
+  fmap not $ readRef (gateRef gate)
+
+-- | Await the gate to be opened if required. If the gate is already open
+-- then the computation returns immediately.
+awaitGateOpened :: Gate -> Process ()
+awaitGateOpened gate =
+  do f <- liftEvent $ readRef (gateRef gate)
+     unless f $
+       do processAwait $ refChanged_ (gateRef gate)
+          awaitGateOpened gate
+
+-- | Await the gate to be closed if required. If the gate is already closed
+-- then the computation returns immediately.
+awaitGateClosed :: Gate -> Process ()
+awaitGateClosed gate =
+  do f <- liftEvent $ readRef (gateRef gate)
+     when f $
+       do processAwait $ refChanged_ (gateRef gate)
+          awaitGateClosed gate
+
+-- | Signal triggered when the state of the gate changes.
+gateChanged_ :: Gate -> Signal ()
+gateChanged_ gate =
+  refChanged_ (gateRef gate)
diff --git a/Simulation/Aivika/Generator.hs b/Simulation/Aivika/Generator.hs
--- a/Simulation/Aivika/Generator.hs
+++ b/Simulation/Aivika/Generator.hs
@@ -1,4 +1,6 @@
 
+{-# LANGUAGE RankNTypes #-}
+
 -- |
 -- Module     : Simulation.Aivika.Generator
 -- Copyright  : Copyright (c) 2009-2015, David Sorokin <david.sorokin@gmail.com>
@@ -12,12 +14,16 @@
 module Simulation.Aivika.Generator 
        (Generator(..),
         GeneratorType(..),
+        DiscretePDF(..),
         newGenerator,
         newRandomGenerator) where
 
 import System.Random
 import Data.IORef
 
+-- | A discrete probability density function.
+type DiscretePDF a = [(a, Double)]
+
 -- | Defines a random number generator.
 data Generator =
   Generator { generateUniform :: Double -> Double -> IO Double,
@@ -26,9 +32,15 @@
               generateUniformInt :: Int -> Int -> IO Int,
               -- ^ Generate an uniform integer random number
               -- with the specified minimum and maximum.
+              generateTriangular :: Double -> Double -> Double -> IO Double,
+              -- ^ Generate a triangular random number
+              -- by the specified minimum, median and maximum.
               generateNormal :: Double -> Double -> IO Double,
               -- ^ Generate the normal random number
               -- with the specified mean and deviation.
+              generateLogNormal :: Double -> Double -> IO Double,
+              -- ^ Generate a random number from the lognormal distribution derived
+              -- from a normal distribution with the specified mean and deviation.
               generateExponential :: Double -> IO Double,
               -- ^ Generate the random number distributed exponentially
               -- with the specified mean (the reciprocal of the rate).
@@ -38,9 +50,28 @@
               generatePoisson :: Double -> IO Int,
               -- ^ Generate the Poisson random number
               -- with the specified mean.
-              generateBinomial :: Double -> Int -> IO Int
+              generateBinomial :: Double -> Int -> IO Int,
               -- ^ Generate the binomial random number
               -- with the specified probability and number of trials.
+              generateGamma :: Double -> Double -> IO Double,
+              -- ^ Generate a random number from the Gamma distribution with
+              -- the specified shape (kappa) and scale (theta, a reciprocal of the rate).
+              --
+              -- The probability density for the Gamma distribution is
+              --
+              -- @f x = x ** (kappa - 1) * exp (- x \/ theta) \/ theta ** kappa * Gamma kappa@
+              generateBeta :: Double -> Double -> IO Double,
+              -- ^ Generate a random number from the Beta distribution by
+              -- the specified shape parameters (alpha and beta).
+              --
+              -- The probability density for the Beta distribution is
+              --
+              -- @f x = x ** (alpha - 1) * (1 - x) ** (beta - 1) \/ B alpha beta@
+              generateWeibull :: Double -> Double -> IO Double,
+              -- ^ Generate a random number from the Weibull distribution by
+              -- the specified shape and scale.
+              generateDiscrete :: forall a. DiscretePDF a -> IO a
+              -- ^ Generate a random value from the specified discrete distribution.
             }
 
 -- | Generate the uniform random number with the specified minimum and maximum.
@@ -69,6 +100,22 @@
          max' = fromIntegral max
      return $ round (min' + x * (max' - min'))
 
+-- | Generate the triangular random number by the specified minimum, median and maximum.
+generateTriangular01 :: IO Double
+                        -- ^ the generator
+                        -> Double
+                        -- ^ minimum
+                        -> Double
+                        -- ^ median
+                        -> Double
+                        -- ^ maximum
+                        -> IO Double
+generateTriangular01 g min median max =
+  do x <- g
+     if x <= (median - min) / (max - min)
+       then return $ min + sqrt ((median - min) * (max - min) * x)
+       else return $ max - sqrt ((max - median) * (max - min) * (1 - x))
+
 -- | Create a normal random number generator with mean 0 and variance 1
 -- by the specified generator of uniform random numbers from 0 to 1.
 newNormalGenerator01 :: IO Double
@@ -169,6 +216,79 @@
                        then loop (n - 1) (acc + 1)
                        else loop (n - 1) acc
 
+-- | Generate a random number from the Gamma distribution using Marsaglia and Tsang method.
+generateGamma01 :: IO Double
+                   -- ^ a normal random number ~ N (0,1)
+                   -> IO Double
+                   -- ^ an uniform random number ~ U (0, 1)
+                   -> Double
+                   -- ^ the shape parameter (kappa) 
+                   -> Double
+                   -- ^ the scale parameter (theta)
+                   -> IO Double
+generateGamma01 gn gu kappa theta
+  | kappa <= 0 = error "The shape parameter (kappa) must be positive: generateGamma01"
+  | kappa > 1  =
+    let d = kappa - 1 / 3
+        c = 1 / sqrt (9 * d)
+        loop =
+          do z <- gn
+             if z <= - (1 / c)
+               then loop
+               else do let v = (1 + c * z) ** 3
+                       u <- gu
+                       if log u > 0.5 * z * z + d - d * v + d * log v
+                         then loop
+                         else return $ d * v * theta
+    in loop
+  | otherwise  =
+    do x <- generateGamma01 gn gu (1 + kappa) theta
+       u <- gu
+       return $ x * u ** (1 / kappa)
+
+-- | Generate a random number from the Beta distribution.
+generateBeta01 :: IO Double
+                  -- ^ a normal random number ~ N (0, 1)
+                  -> IO Double
+                  -- ^ an uniform random number ~ U (0, 1)
+                  -> Double
+                  -- ^ the shape parameter alpha
+                  -> Double
+                  -- ^ the shape parameter beta
+                  -> IO Double
+generateBeta01 gn gu alpha beta =
+  do g1 <- generateGamma01 gn gu alpha 1
+     g2 <- generateGamma01 gn gu beta 1
+     return $ g1 / (g1 + g2)
+
+-- | Generate a random number from the Weibull distribution.
+generateWeibull01 :: IO Double
+                     -- ^ an uniform random number ~ U (0, 1)
+                     -> Double
+                     -- ^ shape
+                     -> Double
+                     -- ^ scale
+                     -> IO Double
+generateWeibull01 g alpha beta =
+  do x <- g
+     return $ beta * (- log x) ** (1 / alpha)
+
+-- | Generate a random value from the specified discrete distribution.
+generateDiscrete01 :: IO Double
+                      -- ^ an uniform random number ~ U (0, 1)
+                      -> DiscretePDF a
+                      -- ^ a discrete probability density function
+                      -> IO a
+generateDiscrete01 g []   = error "Empty PDF: generateDiscrete01"
+generateDiscrete01 g dpdf =
+  do x <- g
+     let loop acc [(a, p)] = a
+         loop acc ((a, p) : dpdf) =
+           if x <= acc + p
+           then a
+           else loop (acc + p) dpdf
+     return $ loop 0 dpdf
+
 -- | Defines a type of the random number generator.
 data GeneratorType = SimpleGenerator
                      -- ^ The simple random number generator.
@@ -208,13 +328,18 @@
 newRandomGenerator01 g =
   do let g1 = g
      g2 <- newNormalGenerator01 g1
-     let g3 mu nu =
-           do x <- g2
-              return $ mu + nu * x
+     let g3 mu nu = do { x <- g2; return $ mu + nu * x }
+         g4 mu nu = do { x <- g2; return $ exp (mu + nu * x) }
      return Generator { generateUniform = generateUniform01 g1,
                         generateUniformInt = generateUniformInt01 g1,
+                        generateTriangular = generateTriangular01 g1,
                         generateNormal = g3,
+                        generateLogNormal = g4,
                         generateExponential = generateExponential01 g1,
                         generateErlang = generateErlang01 g1,
                         generatePoisson = generatePoisson01 g1,
-                        generateBinomial = generateBinomial01 g1 }
+                        generateBinomial = generateBinomial01 g1,
+                        generateGamma = generateGamma01 g2 g1,
+                        generateBeta = generateBeta01 g2 g1,
+                        generateWeibull = generateWeibull01 g1,
+                        generateDiscrete = generateDiscrete01 g1 }
diff --git a/Simulation/Aivika/Net/Random.hs b/Simulation/Aivika/Net/Random.hs
--- a/Simulation/Aivika/Net/Random.hs
+++ b/Simulation/Aivika/Net/Random.hs
@@ -15,12 +15,19 @@
 module Simulation.Aivika.Net.Random
        (randomUniformNet,
         randomUniformIntNet,
+        randomTriangularNet,
         randomNormalNet,
+        randomLogNormalNet,
         randomExponentialNet,
         randomErlangNet,
         randomPoissonNet,
-        randomBinomialNet) where
+        randomBinomialNet,
+        randomGammaNet,
+        randomBetaNet,
+        randomWeibullNet,
+        randomDiscreteNet) where
 
+import Simulation.Aivika.Generator
 import Simulation.Aivika.Process
 import Simulation.Aivika.Process.Random
 import Simulation.Aivika.Net
@@ -48,6 +55,19 @@
   randomUniformIntProcess_ min max
 
 -- | When processing every input element, hold the process
+-- for a random time interval having the triangular distribution.
+randomTriangularNet :: Double
+                       -- ^ the minimum time interval
+                       -> Double
+                       -- ^ the median of the time interval
+                       -> Double
+                       -- ^ the maximum time interval
+                       -> Net a a
+randomTriangularNet min median max =
+  withinNet $
+  randomTriangularProcess_ min median max
+
+-- | When processing every input element, hold the process
 -- for a random time interval distributed normally.
 randomNormalNet :: Double
                    -- ^ the mean time interval
@@ -59,6 +79,19 @@
   randomNormalProcess_ mu nu
          
 -- | When processing every input element, hold the process
+-- for a random time interval having the lognormal distribution.
+randomLogNormalNet :: Double
+                      -- ^ the mean of a normal distribution which
+                      -- this distribution is derived from
+                      -> Double
+                      -- ^ the deviation of a normal distribution which
+                      -- this distribution is derived from
+                      -> Net a a
+randomLogNormalNet mu nu =
+  withinNet $
+  randomLogNormalProcess_ mu nu
+
+-- | When processing every input element, hold the process
 -- for a random time interval distributed exponentially
 -- with the specified mean (the reciprocal of the rate).
 randomExponentialNet :: Double
@@ -101,3 +134,48 @@
 randomBinomialNet prob trials =
   withinNet $
   randomBinomialProcess_ prob trials
+
+-- | When processing every input element, hold the process
+-- for a random time interval having the Gamma distribution
+-- with the specified shape and scale.
+randomGammaNet :: Double
+                  -- ^ the shape
+                  -> Double
+                  -- ^ the scale (a reciprocal of the rate)
+                  -> Net a a
+randomGammaNet kappa theta =
+  withinNet $
+  randomGammaProcess_ kappa theta
+
+-- | When processing every input element, hold the process
+-- for a random time interval having the Beta distribution
+-- with the specified shape parameters (alpha and beta).
+randomBetaNet :: Double
+                 -- ^ shape (alpha)
+                 -> Double
+                 -- ^ shape (beta)
+                 -> Net a a
+randomBetaNet alpha beta =
+  withinNet $
+  randomBetaProcess_ alpha beta
+
+-- | When processing every input element, hold the process
+-- for a random time interval having the Weibull distribution
+-- with the specified shape and scale.
+randomWeibullNet :: Double
+                    -- ^ shape
+                    -> Double
+                    -- ^ scale
+                    -> Net a a
+randomWeibullNet alpha beta =
+  withinNet $
+  randomWeibullProcess_ alpha beta
+
+-- | When processing every input element, hold the process
+-- for a random time interval having the specified discrete distribution.
+randomDiscreteNet :: DiscretePDF Double
+                     -- ^ the discrete probability density function
+                     -> Net a a
+randomDiscreteNet dpdf =
+  withinNet $
+  randomDiscreteProcess_ dpdf
diff --git a/Simulation/Aivika/Parameter/Random.hs b/Simulation/Aivika/Parameter/Random.hs
--- a/Simulation/Aivika/Parameter/Random.hs
+++ b/Simulation/Aivika/Parameter/Random.hs
@@ -22,11 +22,17 @@
 module Simulation.Aivika.Parameter.Random
        (randomUniform,
         randomUniformInt,
+        randomTriangular,
         randomNormal,
+        randomLogNormal,
         randomExponential,
         randomErlang,
         randomPoisson,
         randomBinomial,
+        randomGamma,
+        randomBeta,
+        randomWeibull,
+        randomDiscrete,
         randomTrue,
         randomFalse) where
 
@@ -58,6 +64,16 @@
   let g = runGenerator r
   in generateUniformInt g min max
 
+-- | Computation that generates a new random number from the triangular distribution.
+randomTriangular :: Double     -- ^ minimum
+                    -> Double  -- ^ median
+                    -> Double  -- ^ maximum
+                    -> Parameter Double
+randomTriangular min median max =
+  Parameter $ \r ->
+  let g = runGenerator r
+  in generateTriangular g min median max
+
 -- | Computation that generates a new random number distributed normally.
 randomNormal :: Double     -- ^ mean
                 -> Double  -- ^ deviation
@@ -67,6 +83,19 @@
   let g = runGenerator r
   in generateNormal g mu nu
 
+-- | Computation that generates a new random number from the lognormal distribution.
+randomLogNormal :: Double
+                   -- ^ the mean of a normal distribution
+                   -- which this distribution is derived from
+                   -> Double
+                   -- ^ the deviation of a normal distribution
+                   -- which this distribution is derived from
+                   -> Parameter Double
+randomLogNormal mu nu =
+  Parameter $ \r ->
+  let g = runGenerator r
+  in generateLogNormal g mu nu
+
 -- | Computation that returns a new exponential random number with the specified mean
 -- (the reciprocal of the rate).
 randomExponential :: Double
@@ -121,3 +150,38 @@
 randomFalse p =
   do x <- randomUniform 0 1
      return (x > p)     
+
+-- | Computation that returns a new random number from the Gamma distribution.
+randomGamma :: Double     -- ^ the shape
+               -> Double  -- ^ the scale (a reciprocal of the rate)
+               -> Parameter Double
+randomGamma kappa theta =
+  Parameter $ \r ->
+  let g = runGenerator r
+  in generateGamma g kappa theta
+
+-- | Computation that returns a new random number from the Beta distribution.
+randomBeta :: Double     -- ^ the shape (alpha)
+              -> Double  -- ^ the shape (beta)
+              -> Parameter Double
+randomBeta alpha beta =
+  Parameter $ \r ->
+  let g = runGenerator r
+  in generateBeta g alpha beta
+
+-- | Computation that returns a new random number from the Weibull distribution.
+randomWeibull :: Double     -- ^ shape
+                 -> Double  -- ^ scale
+                 -> Parameter Double
+randomWeibull alpha beta =
+  Parameter $ \r ->
+  let g = runGenerator r
+  in generateWeibull g alpha beta
+
+-- | Computation that returns a new random value from the specified discrete distribution.
+randomDiscrete :: DiscretePDF a -> Parameter a
+randomDiscrete dpdf =
+  Parameter $ \r ->
+  let g = runGenerator r
+  in generateDiscrete g dpdf
+
diff --git a/Simulation/Aivika/PriorityQueue.hs b/Simulation/Aivika/PriorityQueue.hs
--- a/Simulation/Aivika/PriorityQueue.hs
+++ b/Simulation/Aivika/PriorityQueue.hs
@@ -19,6 +19,8 @@
         enqueue, 
         dequeue, 
         queueFront,
+        queueDelete,
+        queueDeleteBy,
         remove,
         removeBy) where 
 
@@ -26,6 +28,8 @@
 import Data.Array.MArray.Safe
 import Data.Array.IO.Safe
 import Data.IORef
+import Data.Maybe
+
 import Control.Monad
 
 -- | The 'PriorityQueue' type represents an imperative heap-based 
@@ -171,25 +175,26 @@
 -- indicating whether the element was actually removed.
 --
 -- Note that unlike other functions it has complexity O(n).
-remove :: Eq a => PriorityQueue a -> a -> IO Bool
-remove pq a = removeBy pq (== a)
+queueDelete :: Eq a => PriorityQueue a -> a -> IO Bool
+queueDelete pq a = fmap isJust $ queueDeleteBy pq (== a)
 
--- | Remove an element satisfying the predicate and return a computation of the flag
--- indicating whether the element was actually removed.
+-- | Remove an element satisfying the predicate and return a computation of
+-- the element if found.
 --
 -- Note that unlike other functions it has complexity O(n).
-removeBy :: PriorityQueue a -> (a -> Bool) -> IO Bool
-removeBy pq pred =
-  do index <- indexBy pq pred
+queueDeleteBy :: PriorityQueue a -> (a -> Bool) -> IO (Maybe a)
+queueDeleteBy pq pred =
+  do index <- queueIndexBy pq pred
      if index < 0
-       then return False
+       then return Nothing
        else do size <- readIORef (pqSize pq)
                when (size == 0) $
-                 error "Internal error in the priority queue implementation: removeBy"
+                 error "Internal error in the priority queue implementation: queueDeleteBy"
                let i = size - 1
                writeIORef (pqSize pq) i
                keys <- readIORef (pqKeys pq)
                vals <- readIORef (pqVals pq)
+               x <- readArray vals index
                k <- readArray keys i
                v <- readArray vals i
                let k0 = 0.0
@@ -198,11 +203,11 @@
                writeArray vals i v0
                when (i > 0) $
                  siftDown keys vals i index k v
-               return True
+               return (Just x)
      
 -- | Return the index of the item satisfying the predicate or -1.     
-indexBy :: PriorityQueue a -> (a -> Bool) -> IO Int
-indexBy pq pred =
+queueIndexBy :: PriorityQueue a -> (a -> Bool) -> IO Int
+queueIndexBy pq pred =
   do size <- readIORef (pqSize pq)
      vals <- readIORef (pqVals pq)
      let loop index =
@@ -213,3 +218,13 @@
                      then return index
                      else loop $ index + 1
      loop 0
+
+-- | Use 'queueDelete' instead.
+remove :: Eq a => PriorityQueue a -> a -> IO Bool
+{-# DEPRECATED remove "Use queueDelete instead." #-}
+remove = queueDelete
+
+-- | Use 'queueDeleteBy' instead.
+removeBy :: PriorityQueue a -> (a -> Bool) -> IO Bool
+{-# DEPRECATED removeBy "Use queueDeleteBy instead." #-}
+removeBy pq pred = fmap isJust $ queueDeleteBy pq pred
diff --git a/Simulation/Aivika/Process/Random.hs b/Simulation/Aivika/Process/Random.hs
--- a/Simulation/Aivika/Process/Random.hs
+++ b/Simulation/Aivika/Process/Random.hs
@@ -17,8 +17,12 @@
         randomUniformProcess_,
         randomUniformIntProcess,
         randomUniformIntProcess_,
+        randomTriangularProcess,
+        randomTriangularProcess_,
         randomNormalProcess,
         randomNormalProcess_,
+        randomLogNormalProcess,
+        randomLogNormalProcess_,
         randomExponentialProcess,
         randomExponentialProcess_,
         randomErlangProcess,
@@ -26,11 +30,20 @@
         randomPoissonProcess,
         randomPoissonProcess_,
         randomBinomialProcess,
-        randomBinomialProcess_) where
+        randomBinomialProcess_,
+        randomGammaProcess,
+        randomGammaProcess_,
+        randomBetaProcess,
+        randomBetaProcess_,
+        randomWeibullProcess,
+        randomWeibullProcess_,
+        randomDiscreteProcess,
+        randomDiscreteProcess_) where
 
 import Control.Monad
 import Control.Monad.Trans
 
+import Simulation.Aivika.Generator
 import Simulation.Aivika.Parameter
 import Simulation.Aivika.Parameter.Random
 import Simulation.Aivika.Process
@@ -81,6 +94,33 @@
   do t <- liftParameter $ randomUniformInt min max
      holdProcess $ fromIntegral t
 
+-- | Hold the process for a random time interval having the triangular distribution.
+randomTriangularProcess :: Double
+                           -- ^ the minimum time interval
+                           -> Double
+                           -- ^ a median of the time interval
+                           -> Double
+                           -- ^ the maximum time interval
+                           -> Process Double
+                           -- ^ a computation of the time interval
+                           -- for which the process was actually held
+randomTriangularProcess min median max =
+  do t <- liftParameter $ randomTriangular min median max
+     holdProcess t
+     return t
+
+-- | Hold the process for a random time interval having the triangular distribution.
+randomTriangularProcess_ :: Double
+                            -- ^ the minimum time interval
+                            -> Double
+                            -- ^ a median of the time interval
+                            -> Double
+                            -- ^ the maximum time interval
+                            -> Process ()
+randomTriangularProcess_ min median max =
+  do t <- liftParameter $ randomTriangular min median max
+     holdProcess t
+
 -- | Hold the process for a random time interval distributed normally.
 randomNormalProcess :: Double
                        -- ^ the mean time interval
@@ -105,7 +145,34 @@
   do t <- liftParameter $ randomNormal mu nu
      when (t > 0) $
        holdProcess t
-         
+
+-- | Hold the process for a random time interval having the lognormal distribution.
+randomLogNormalProcess :: Double
+                          -- ^ the mean for a normal distribution
+                          -- which this distribution is derived from
+                          -> Double
+                          -- ^ the deviation for a normal distribution
+                          -- which this distribution is derived from
+                          -> Process Double
+                          -- ^ a computation of the time interval
+                          -- for which the process was actually held
+randomLogNormalProcess mu nu =
+  do t <- liftParameter $ randomLogNormal mu nu
+     holdProcess t
+     return t
+
+-- | Hold the process for a random time interval having the lognormal distribution.
+randomLogNormalProcess_ :: Double
+                           -- ^ the mean for a normal distribution
+                           -- which this distribution is derived from
+                           -> Double
+                           -- ^ the deviation for a normal distribution
+                           -- which this distribution is derived from
+                           -> Process ()
+randomLogNormalProcess_ mu nu =
+  do t <- liftParameter $ randomLogNormal mu nu
+     holdProcess t
+
 -- | Hold the process for a random time interval distributed exponentially
 -- with the specified mean (the reciprocal of the rate).
 randomExponentialProcess :: Double
@@ -197,3 +264,97 @@
 randomBinomialProcess_ prob trials =
   do t <- liftParameter $ randomBinomial prob trials
      holdProcess $ fromIntegral t
+
+-- | Hold the process for a random time interval having the Gamma distribution
+-- with the specified shape and scale.
+randomGammaProcess :: Double
+                      -- ^ the shape
+                      -> Double
+                      -- ^ the scale (a reciprocal of the rate)
+                      -> Process Double
+                      -- ^ a computation of the time interval
+                      -- for which the process was actually held
+randomGammaProcess kappa theta =
+  do t <- liftParameter $ randomGamma kappa theta
+     holdProcess t
+     return t
+
+-- | Hold the process for a random time interval having the Gamma distribution
+-- with the specified shape and scale.
+randomGammaProcess_ :: Double
+                       -- ^ the shape
+                       -> Double
+                       -- ^ the scale (a reciprocal of the rate)
+                       -> Process ()
+randomGammaProcess_ kappa theta =
+  do t <- liftParameter $ randomGamma kappa theta
+     holdProcess t
+
+-- | Hold the process for a random time interval having the Beta distribution
+-- with the specified shape parameters (alpha and beta).
+randomBetaProcess :: Double
+                     -- ^ the shape (alpha)
+                     -> Double
+                     -- ^ the shape (beta)
+                     -> Process Double
+                     -- ^ a computation of the time interval
+                     -- for which the process was actually held
+randomBetaProcess alpha beta =
+  do t <- liftParameter $ randomBeta alpha beta
+     holdProcess t
+     return t
+
+-- | Hold the process for a random time interval having the Beta distribution
+-- with the specified shape parameters (alpha and beta).
+randomBetaProcess_ :: Double
+                      -- ^ the shape (alpha)
+                      -> Double
+                      -- ^ the shape (beta)
+                      -> Process ()
+randomBetaProcess_ alpha beta =
+  do t <- liftParameter $ randomBeta alpha beta
+     holdProcess t
+
+-- | Hold the process for a random time interval having the Weibull distribution
+-- with the specified shape and scale.
+randomWeibullProcess :: Double
+                        -- ^ the shape
+                        -> Double
+                        -- ^ the scale
+                        -> Process Double
+                        -- ^ a computation of the time interval
+                        -- for which the process was actually held
+randomWeibullProcess alpha beta =
+  do t <- liftParameter $ randomWeibull alpha beta
+     holdProcess t
+     return t
+
+-- | Hold the process for a random time interval having the Weibull distribution
+-- with the specified shape and scale.
+randomWeibullProcess_ :: Double
+                         -- ^ the shape
+                         -> Double
+                         -- ^ the scale
+                         -> Process ()
+randomWeibullProcess_ alpha beta =
+  do t <- liftParameter $ randomWeibull alpha beta
+     holdProcess t
+
+-- | Hold the process for a random time interval having the specified discrete distribution.
+randomDiscreteProcess :: DiscretePDF Double
+                         -- ^ the discrete probability density function
+                         -> Process Double
+                         -- ^ a computation of the time interval
+                         -- for which the process was actually held
+randomDiscreteProcess dpdf =
+  do t <- liftParameter $ randomDiscrete dpdf
+     holdProcess t
+     return t
+
+-- | Hold the process for a random time interval having the specified discrete distribution.
+randomDiscreteProcess_ :: DiscretePDF Double
+                          -- ^ the discrete probability density function
+                          -> Process ()
+randomDiscreteProcess_ dpdf =
+  do t <- liftParameter $ randomDiscrete dpdf
+     holdProcess t
diff --git a/Simulation/Aivika/Processor/Random.hs b/Simulation/Aivika/Processor/Random.hs
--- a/Simulation/Aivika/Processor/Random.hs
+++ b/Simulation/Aivika/Processor/Random.hs
@@ -15,12 +15,19 @@
 module Simulation.Aivika.Processor.Random
        (randomUniformProcessor,
         randomUniformIntProcessor,
+        randomTriangularProcessor,
         randomNormalProcessor,
+        randomLogNormalProcessor,
         randomExponentialProcessor,
         randomErlangProcessor,
         randomPoissonProcessor,
-        randomBinomialProcessor) where
+        randomBinomialProcessor,
+        randomGammaProcessor,
+        randomBetaProcessor,
+        randomWeibullProcessor,
+        randomDiscreteProcessor) where
 
+import Simulation.Aivika.Generator
 import Simulation.Aivika.Process
 import Simulation.Aivika.Process.Random
 import Simulation.Aivika.Processor
@@ -48,6 +55,19 @@
   randomUniformIntProcess_ min max
 
 -- | When processing every input element, hold the process
+-- for a random time interval having the triangular distribution.
+randomTriangularProcessor :: Double
+                             -- ^ the minimum time interval
+                             -> Double
+                             -- ^ the median of the time interval
+                             -> Double
+                             -- ^ the maximum time interval
+                             -> Processor a a
+randomTriangularProcessor min median max =
+  withinProcessor $
+  randomTriangularProcess_ min median max
+
+-- | When processing every input element, hold the process
 -- for a random time interval distributed normally.
 randomNormalProcessor :: Double
                          -- ^ the mean time interval
@@ -59,6 +79,19 @@
   randomNormalProcess_ mu nu
          
 -- | When processing every input element, hold the process
+-- for a random time interval having the lognormal distribution.
+randomLogNormalProcessor :: Double
+                            -- ^ the mean for a normal distribution
+                            -- which this distribution is derived from
+                            -> Double
+                            -- ^ the deviation for a normal distribution
+                            -- which this distribution is derived from
+                            -> Processor a a
+randomLogNormalProcessor mu nu =
+  withinProcessor $
+  randomLogNormalProcess_ mu nu
+
+-- | When processing every input element, hold the process
 -- for a random time interval distributed exponentially
 -- with the specified mean (the reciprocal of the rate).
 randomExponentialProcessor :: Double
@@ -101,3 +134,48 @@
 randomBinomialProcessor prob trials =
   withinProcessor $
   randomBinomialProcess_ prob trials
+
+-- | When processing every input element, hold the process
+-- for a random time interval having the Gamma distribution
+-- with the specified shape and scale.
+randomGammaProcessor :: Double
+                        -- ^ the shape
+                        -> Double
+                        -- ^ the scale (a reciprocal of the rate)
+                        -> Processor a a
+randomGammaProcessor kappa theta =
+  withinProcessor $
+  randomGammaProcess_ kappa theta
+
+-- | When processing every input element, hold the process
+-- for a random time interval having the Beta distribution
+-- with the specified shape parameters (alpha and beta).
+randomBetaProcessor :: Double
+                       -- ^ shape (alpha)
+                       -> Double
+                       -- ^ shape (beta)
+                       -> Processor a a
+randomBetaProcessor alpha beta =
+  withinProcessor $
+  randomBetaProcess_ alpha beta
+
+-- | When processing every input element, hold the process
+-- for a random time interval having the Weibull distribution
+-- with the specified shape and scale.
+randomWeibullProcessor :: Double
+                          -- ^ shape
+                          -> Double
+                          -- ^ scale
+                          -> Processor a a
+randomWeibullProcessor alpha beta =
+  withinProcessor $
+  randomWeibullProcess_ alpha beta
+
+-- | When processing every input element, hold the process
+-- for a random time interval having the specified discrete distribution.
+randomDiscreteProcessor :: DiscretePDF Double
+                           -- ^ the discrete probability density function
+                           -> Processor a a
+randomDiscreteProcessor dpdf =
+  withinProcessor $
+  randomDiscreteProcess_ dpdf
diff --git a/Simulation/Aivika/Queue.hs b/Simulation/Aivika/Queue.hs
--- a/Simulation/Aivika/Queue.hs
+++ b/Simulation/Aivika/Queue.hs
@@ -63,6 +63,11 @@
         enqueueOrLost_,
         enqueueWithStoringPriorityOrLost,
         enqueueWithStoringPriorityOrLost_,
+        queueDelete,
+        queueDelete_,
+        queueDeleteBy,
+        queueDeleteBy_,
+        clearQueue,
         -- * Awaiting
         waitWhileFullQueue,
         -- * Summary
@@ -107,6 +112,7 @@
 
 import Data.IORef
 import Data.Monoid
+import Data.Maybe
 
 import Control.Monad
 import Control.Monad.Trans
@@ -487,7 +493,7 @@
 -- | Return the total wait time from the time at which the enqueueing operation
 -- was initiated to the time at which the item was dequeued.
 --
--- In some sense, @queueTotalWaitTime == queueInputWaitTime + queueWaitTime@.
+-- In some sense, @queueTotalWaitTime == enqueueWaitTime + queueWaitTime@.
 --
 -- See also 'queueTotalWaitTimeChanged' and 'queueTotalWaitTimeChanged_'.
 queueTotalWaitTime :: Queue si sm so a -> Event (SamplingStats Double)
@@ -607,6 +613,77 @@
                fmap Just $ dequeueExtract q t
        else return Nothing
 
+-- | Remove the item from the queue and return a flag indicating
+-- whether the item was found and actually removed.
+queueDelete :: (Eq a,
+                DequeueStrategy si,
+                DeletingQueueStrategy sm,
+                DequeueStrategy so)
+               => Queue si sm so a
+               -- ^ the queue
+               -> a
+               -- ^ the item to remove from the queue
+               -> Event Bool
+               -- ^ whether the item was found and removed
+queueDelete q a = fmap isJust $ queueDeleteBy q (== a)
+
+-- | Remove the specified item from the queue.
+queueDelete_ :: (Eq a,
+                 DequeueStrategy si,
+                 DeletingQueueStrategy sm,
+                 DequeueStrategy so)
+                => Queue si sm so a
+                -- ^ the queue
+                -> a
+                -- ^ the item to remove from the queue
+                -> Event ()
+queueDelete_ q a = fmap (const ()) $ queueDeleteBy q (== a)
+
+-- | Remove an item satisfying the specified predicate and return the item if found.
+queueDeleteBy :: (DequeueStrategy si,
+                  DeletingQueueStrategy sm,
+                  DequeueStrategy so)
+                 => Queue si sm so a
+                 -- ^ the queue
+                 -> (a -> Bool)
+                 -- ^ the predicate
+                 -> Event (Maybe a)
+queueDeleteBy q pred =
+  do x <- tryRequestResourceWithinEvent (dequeueRes q)
+     if x
+       then do i <- strategyQueueDeleteBy (queueStore q) (pred . itemValue)
+               case i of
+                 Nothing ->
+                   do releaseResourceWithinEvent (dequeueRes q)
+                      return Nothing
+                 Just i ->
+                   do t <- dequeueRequest q
+                      fmap Just $ dequeuePostExtract q t i
+       else return Nothing
+               
+-- | Remove an item satisfying the specified predicate.
+queueDeleteBy_ :: (DequeueStrategy si,
+                   DeletingQueueStrategy sm,
+                   DequeueStrategy so)
+                  => Queue si sm so a
+                  -- ^ the queue
+                  -> (a -> Bool)
+                  -- ^ the predicate
+                  -> Event ()
+queueDeleteBy_ q pred = fmap (const ()) $ queueDeleteBy q pred
+
+-- | Clear the queue immediately.
+clearQueue :: (DequeueStrategy si,
+               DequeueStrategy sm)
+              => Queue si sm so a
+              -- ^ the queue
+              -> Event ()
+clearQueue q =
+  do x <- tryDequeue q
+     case x of
+       Nothing -> return ()
+       Just a  -> clearQueue q
+              
 -- | Enqueue the item suspending the process if the queue is full.  
 enqueue :: (EnqueueStrategy si,
             EnqueueStrategy sm,
@@ -921,7 +998,23 @@
   Event $ \p ->
   do i <- invokeEvent p $
           strategyDequeue (queueStore q)
-     c <- readIORef (queueCountRef q)
+     invokeEvent p $
+       dequeuePostExtract q t' i
+
+-- | A post action after extracting the item by the dequeuing request.  
+dequeuePostExtract :: (DequeueStrategy si,
+                       DequeueStrategy sm)
+                      => Queue si sm so a
+                      -- ^ the queue
+                      -> Double
+                      -- ^ the time of the dequeuing request
+                      -> QueueItem a
+                      -- ^ the item to dequeue
+                      -> Event a
+                      -- ^ the dequeued value
+dequeuePostExtract q t' i =
+  Event $ \p ->
+  do c <- readIORef (queueCountRef q)
      let c' = c - 1
          t  = pointTime p
      c' `seq` writeIORef (queueCountRef q) c'
diff --git a/Simulation/Aivika/Queue/Infinite.hs b/Simulation/Aivika/Queue/Infinite.hs
--- a/Simulation/Aivika/Queue/Infinite.hs
+++ b/Simulation/Aivika/Queue/Infinite.hs
@@ -43,6 +43,11 @@
         tryDequeue,
         enqueue,
         enqueueWithStoringPriority,
+        queueDelete,
+        queueDelete_,
+        queueDeleteBy,
+        queueDeleteBy_,
+        clearQueue,
         -- * Summary
         queueSummary,
         -- * Derived Signals for Properties
@@ -71,6 +76,7 @@
 
 import Data.IORef
 import Data.Monoid
+import Data.Maybe
 
 import Control.Monad
 import Control.Monad.Trans
@@ -408,6 +414,72 @@
                fmap Just $ dequeueExtract q t
        else return Nothing
 
+-- | Remove the item from the queue and return a flag indicating
+-- whether the item was found and actually removed.
+queueDelete :: (Eq a,
+                DeletingQueueStrategy sm,
+                DequeueStrategy so)
+               => Queue sm so a
+               -- ^ the queue
+               -> a
+               -- ^ the item to remove from the queue
+               -> Event Bool
+               -- ^ whether the item was found and removed
+queueDelete q a = fmap isJust $ queueDeleteBy q (== a)
+
+-- | Remove the specified item from the queue.
+queueDelete_ :: (Eq a,
+                 DeletingQueueStrategy sm,
+                 DequeueStrategy so)
+                => Queue sm so a
+                -- ^ the queue
+                -> a
+                -- ^ the item to remove from the queue
+                -> Event ()
+queueDelete_ q a = fmap (const ()) $ queueDeleteBy q (== a)
+
+-- | Remove an item satisfying the specified predicate and return the item if found.
+queueDeleteBy :: (DeletingQueueStrategy sm,
+                  DequeueStrategy so)
+                 => Queue sm so a
+                 -- ^ the queue
+                 -> (a -> Bool)
+                 -- ^ the predicate
+                 -> Event (Maybe a)
+queueDeleteBy q pred =
+  do x <- tryRequestResourceWithinEvent (dequeueRes q)
+     if x
+       then do i <- strategyQueueDeleteBy (queueStore q) (pred . itemValue)
+               case i of
+                 Nothing ->
+                   do releaseResourceWithinEvent (dequeueRes q)
+                      return Nothing
+                 Just i ->
+                   do t <- dequeueRequest q
+                      fmap Just $ dequeuePostExtract q t i
+       else return Nothing
+               
+-- | Remove an item satisfying the specified predicate.
+queueDeleteBy_ :: (DeletingQueueStrategy sm,
+                   DequeueStrategy so)
+                  => Queue sm so a
+                  -- ^ the queue
+                  -> (a -> Bool)
+                  -- ^ the predicate
+                  -> Event ()
+queueDeleteBy_ q pred = fmap (const ()) $ queueDeleteBy q pred
+
+-- | Clear the queue immediately.
+clearQueue :: DequeueStrategy sm
+              => Queue sm so a
+              -- ^ the queue
+              -> Event ()
+clearQueue q =
+  do x <- tryDequeue q
+     case x of
+       Nothing -> return ()
+       Just a  -> clearQueue q
+
 -- | Enqueue the item.  
 enqueue :: (EnqueueStrategy sm,
             DequeueStrategy so)
@@ -520,7 +592,22 @@
   Event $ \p ->
   do i <- invokeEvent p $
           strategyDequeue (queueStore q)
-     c <- readIORef (queueCountRef q)
+     invokeEvent p $
+       dequeuePostExtract q t' i
+
+-- | A post action after extracting the item by the dequeuing request.  
+dequeuePostExtract :: DequeueStrategy sm
+                      => Queue sm so a
+                      -- ^ the queue
+                      -> Double
+                      -- ^ the time of the dequeuing request
+                      -> QueueItem a
+                      -- ^ the item to dequeue
+                      -> Event a
+                      -- ^ the dequeued value
+dequeuePostExtract q t' i =
+  Event $ \p ->
+  do c <- readIORef (queueCountRef q)
      let c' = c - 1
          t  = pointTime p
      c' `seq` writeIORef (queueCountRef q) c'
diff --git a/Simulation/Aivika/QueueStrategy.hs b/Simulation/Aivika/QueueStrategy.hs
--- a/Simulation/Aivika/QueueStrategy.hs
+++ b/Simulation/Aivika/QueueStrategy.hs
@@ -15,6 +15,7 @@
 
 import System.Random
 import Control.Monad.Trans
+import Data.Maybe
 
 import Simulation.Aivika.Simulation
 import Simulation.Aivika.Event
@@ -73,6 +74,28 @@
                                  -> Event ()
                                  -- ^ the action of enqueuing
 
+-- | Defines a strategy with support of the deleting operation.
+class DequeueStrategy s => DeletingQueueStrategy s where
+
+  -- | Remove the element and return a flag indicating whether
+  -- the element was found and removed.
+  strategyQueueDelete :: Eq i
+                         => StrategyQueue s i
+                         -- ^ the queue
+                         -> i
+                         -- ^ the element
+                         -> Event Bool
+                         -- ^ whether the element was found and removed
+  strategyQueueDelete s i = fmap isJust $ strategyQueueDeleteBy s (== i)
+
+  -- | Remove an element satisfying the predicate and return the element if found.
+  strategyQueueDeleteBy :: StrategyQueue s i
+                           -- ^ the queue
+                           -> (i -> Bool)
+                           -- ^ the predicate
+                           -> Event (Maybe i)
+                           -- ^ the element if it was found and removed
+
 -- | Strategy: First Come - First Served (FCFS).
 data FCFS = FCFS deriving (Eq, Ord, Show)
 
@@ -109,6 +132,11 @@
 
   strategyEnqueue (FCFSQueue q) i = liftIO $ listAddLast q i
 
+-- | An implementation of the 'FCFS' queue strategy.
+instance DeletingQueueStrategy FCFS where
+
+  strategyQueueDeleteBy (FCFSQueue q) p = liftIO $ listRemoveBy q p
+
 -- | An implementation of the 'LCFS' queue strategy.
 instance QueueStrategy LCFS where
 
@@ -133,6 +161,11 @@
 
   strategyEnqueue (LCFSQueue q) i = liftIO $ listInsertFirst q i
 
+-- | An implementation of the 'LCFS' queue strategy.
+instance DeletingQueueStrategy LCFS where
+
+  strategyQueueDeleteBy (LCFSQueue q) p = liftIO $ listRemoveBy q p
+
 -- | An implementation of the 'StaticPriorities' queue strategy.
 instance QueueStrategy StaticPriorities where
 
@@ -157,6 +190,11 @@
 
   strategyEnqueueWithPriority (StaticPriorityQueue q) p i = liftIO $ PQ.enqueue q p i
 
+-- | An implementation of the 'StaticPriorities' queue strategy.
+instance DeletingQueueStrategy StaticPriorities where
+
+  strategyQueueDeleteBy (StaticPriorityQueue q) p = liftIO $ PQ.queueDeleteBy q p
+
 -- | An implementation of the 'SIRO' queue strategy.
 instance QueueStrategy SIRO where
 
@@ -185,3 +223,8 @@
 instance EnqueueStrategy SIRO where
 
   strategyEnqueue (SIROQueue q) i = liftIO $ V.appendVector q i
+
+-- | An implementation of the 'SIRO' queue strategy.
+instance DeletingQueueStrategy SIRO where
+
+  strategyQueueDeleteBy (SIROQueue q) p = liftIO $ V.vectorDeleteBy q p
diff --git a/Simulation/Aivika/Resource/Preemption.hs b/Simulation/Aivika/Resource/Preemption.hs
--- a/Simulation/Aivika/Resource/Preemption.hs
+++ b/Simulation/Aivika/Resource/Preemption.hs
@@ -46,6 +46,7 @@
 
 import Data.IORef
 import Data.Monoid
+import Data.Maybe
 
 import Control.Monad
 import Control.Monad.Trans
@@ -308,7 +309,7 @@
   Process $ \pid ->
   Cont $ \c ->
   Event $ \p ->
-  do f <- PQ.removeBy (resourceActingQueue r) (\item -> actingItemId item == pid)
+  do f <- fmap isJust $ PQ.queueDeleteBy (resourceActingQueue r) (\item -> actingItemId item == pid)
      if f
        then do invokeEvent p $ updateResourceUtilisationCount r (-1)
                invokeEvent p $ releaseResource' r
diff --git a/Simulation/Aivika/Resource/Preemption/Base.hs b/Simulation/Aivika/Resource/Preemption/Base.hs
--- a/Simulation/Aivika/Resource/Preemption/Base.hs
+++ b/Simulation/Aivika/Resource/Preemption/Base.hs
@@ -32,6 +32,7 @@
         alterResourceCount) where
 
 import Data.IORef
+import Data.Maybe
 
 import Control.Monad
 import Control.Monad.Trans
@@ -186,7 +187,7 @@
   Process $ \pid ->
   Cont $ \c ->
   Event $ \p ->
-  do f <- PQ.removeBy (resourceActingQueue r) (\item -> actingItemId item == pid)
+  do f <- fmap isJust $ PQ.queueDeleteBy (resourceActingQueue r) (\item -> actingItemId item == pid)
      if f
        then do invokeEvent p $ releaseResource' r
                invokeEvent p $ resumeCont c ()
diff --git a/Simulation/Aivika/Results/Transform.hs b/Simulation/Aivika/Results/Transform.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Results/Transform.hs
@@ -0,0 +1,610 @@
+
+{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}
+
+-- |
+-- Module     : Simulation.Aivika.Results.Transform
+-- Copyright  : Copyright (c) 2009-2015, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.1
+--
+-- The module defines useful result transformations that can
+-- be used in simulation experiments.
+--
+module Simulation.Aivika.Results.Transform
+       (-- * Basic Class Type
+        ResultTransformer(..),
+        -- * Sampling Statistics
+        SamplingStats(..),
+        samplingStatsCount,
+        samplingStatsMin,
+        samplingStatsMax,
+        samplingStatsMean,
+        samplingStatsMean2,
+        samplingStatsVariance,
+        samplingStatsDeviation,
+        -- * Time-dependent Statistics
+        TimingStats(..),
+        timingStatsCount,
+        timingStatsMin,
+        timingStatsMax,
+        timingStatsMean,
+        timingStatsVariance,
+        timingStatsDeviation,
+        timingStatsMinTime,
+        timingStatsMaxTime,
+        timingStatsStartTime,
+        timingStatsLastTime,
+        timingStatsSum,
+        timingStatsSum2,
+        -- * Sampling-based Counter
+        SamplingCounter(..),
+        samplingCounterValue,
+        samplingCounterStats,
+        -- * Time-dependent Counter
+        TimingCounter(..),
+        timingCounterValue,
+        timingCounterStats,
+        -- * Queue
+        Queue(..),
+        enqueueStrategy,
+        enqueueStoringStrategy,
+        dequeueStrategy,
+        queueNull,
+        queueFull,
+        queueMaxCount,
+        queueCount,
+        queueCountStats,
+        enqueueCount,
+        enqueueLostCount,
+        enqueueStoreCount,
+        dequeueCount,
+        dequeueExtractCount,
+        queueLoadFactor,
+        enqueueRate,
+        enqueueStoreRate,
+        dequeueRate,
+        dequeueExtractRate,
+        queueWaitTime,
+        queueTotalWaitTime,
+        enqueueWaitTime,
+        dequeueWaitTime,
+        queueRate,
+        -- * Arrival Timer
+        ArrivalTimer(..),
+        arrivalProcessingTime,
+        -- * Server
+        Server(..),
+        serverInitState,
+        serverState,
+        serverTotalInputWaitTime,
+        serverTotalProcessingTime,
+        serverTotalOutputWaitTime,
+        serverTotalPreemptionTime,
+        serverInputWaitTime,
+        serverProcessingTime,
+        serverOutputWaitTime,
+        serverPreemptionTime,
+        serverInputWaitFactor,
+        serverProcessingFactor,
+        serverOutputWaitFactor,
+        serverPreemptionFactor,
+        -- * Activity
+        Activity(..),
+        activityInitState,
+        activityState,
+        activityTotalUtilisationTime,
+        activityTotalIdleTime,
+        activityTotalPreemptionTime,
+        activityUtilisationTime,
+        activityIdleTime,
+        activityPreemptionTime,
+        activityUtilisationFactor,
+        activityIdleFactor,
+        activityPreemptionFactor,
+        -- * Resource
+        Resource(..),
+        resourceCount,
+        resourceCountStats,
+        resourceUtilisationCount,
+        resourceUtilisationCountStats,
+        resourceQueueCount,
+        resourceQueueCountStats,
+        resourceTotalWaitTime,
+        resourceWaitTime) where
+
+import Control.Arrow
+
+import Simulation.Aivika.Results
+import Simulation.Aivika.Results.Locale
+
+-- | Something that can transform the results.
+class ResultTransformer a where
+
+  -- | Return the result transform.
+  tr :: a -> ResultTransform
+
+instance ResultTransformer ResultTransform where
+  tr = id
+
+-- | Represents a statistics based upon observations.
+newtype SamplingStats = SamplingStats ResultTransform
+
+instance ResultTransformer SamplingStats where
+  tr (SamplingStats a) = a
+
+-- | The total number of samples.
+samplingStatsCount :: SamplingStats -> ResultTransform
+samplingStatsCount (SamplingStats a) =
+  a >>> expandResults >>> resultById SamplingStatsCountId
+
+-- | The minimum value among the samples.
+samplingStatsMin :: SamplingStats -> ResultTransform
+samplingStatsMin (SamplingStats a) =
+  a >>> expandResults >>> resultById SamplingStatsMinId
+
+-- | The maximum value among the samples.
+samplingStatsMax :: SamplingStats -> ResultTransform
+samplingStatsMax (SamplingStats a) =
+  a >>> expandResults >>> resultById SamplingStatsMaxId
+  
+-- | The average value.
+samplingStatsMean :: SamplingStats -> ResultTransform
+samplingStatsMean (SamplingStats a) =
+  a >>> expandResults >>> resultById SamplingStatsMeanId
+
+-- | The average square value.
+samplingStatsMean2 :: SamplingStats -> ResultTransform
+samplingStatsMean2 (SamplingStats a) =
+  a >>> expandResults >>> resultById SamplingStatsMean2Id
+
+-- | Return tha variance.
+samplingStatsVariance :: SamplingStats -> ResultTransform
+samplingStatsVariance (SamplingStats a) =
+  a >>> expandResults >>> resultById SamplingStatsVarianceId
+
+-- | Return the deviation.
+samplingStatsDeviation :: SamplingStats -> ResultTransform
+samplingStatsDeviation (SamplingStats a) =
+  a >>> expandResults >>> resultById SamplingStatsDeviationId
+
+-- | A counter for which the statistics is collected too.
+newtype SamplingCounter = SamplingCounter ResultTransform
+
+instance ResultTransformer SamplingCounter where
+  tr (SamplingCounter a) = a
+
+-- | The counter value.
+samplingCounterValue :: SamplingCounter -> ResultTransform
+samplingCounterValue (SamplingCounter a) =
+  a >>> resultById SamplingCounterValueId
+
+-- | The counter statistics.
+samplingCounterStats :: SamplingCounter -> SamplingStats
+samplingCounterStats (SamplingCounter a) =
+  SamplingStats (a >>> resultById SamplingCounterStatsId)
+
+-- | The time-dependent statistics.
+newtype TimingStats = TimingStats ResultTransform
+
+instance ResultTransformer TimingStats where
+  tr (TimingStats a) = a
+
+-- | Return the number of samples.
+timingStatsCount :: TimingStats -> ResultTransform
+timingStatsCount (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsCountId
+
+-- | Return the minimum value.
+timingStatsMin :: TimingStats -> ResultTransform
+timingStatsMin (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsMinId
+
+-- | Return the maximum value.
+timingStatsMax :: TimingStats -> ResultTransform
+timingStatsMax (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsMaxId
+
+-- | Return the average value.
+timingStatsMean :: TimingStats -> ResultTransform
+timingStatsMean (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsMeanId
+
+-- | Return the variance.
+timingStatsVariance :: TimingStats -> ResultTransform
+timingStatsVariance (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsVarianceId
+
+-- | Return the deviation.
+timingStatsDeviation :: TimingStats -> ResultTransform
+timingStatsDeviation (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsDeviationId
+
+-- | Return the time at which the minimum is attained.
+timingStatsMinTime :: TimingStats -> ResultTransform
+timingStatsMinTime (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsMinTimeId
+
+-- | Return the time at which the maximum is attained.
+timingStatsMaxTime :: TimingStats -> ResultTransform
+timingStatsMaxTime (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsMaxTimeId
+
+-- | Return the start time of sampling.
+timingStatsStartTime :: TimingStats -> ResultTransform
+timingStatsStartTime (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsStartTimeId
+
+-- | Return the last time of sampling.
+timingStatsLastTime :: TimingStats -> ResultTransform
+timingStatsLastTime (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsLastTimeId
+
+-- | Return the sum of values.
+timingStatsSum :: TimingStats -> ResultTransform
+timingStatsSum (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsSumId
+
+-- | Return the sum of square values.
+timingStatsSum2 :: TimingStats -> ResultTransform
+timingStatsSum2 (TimingStats a) =
+  a >>> expandResults >>> resultById TimingStatsSum2Id
+
+-- | A time-dependent counter that collects the statistics too.
+newtype TimingCounter = TimingCounter ResultTransform
+
+instance ResultTransformer TimingCounter where
+  tr (TimingCounter a) = a
+
+-- | The counter value.
+timingCounterValue :: TimingCounter -> ResultTransform
+timingCounterValue (TimingCounter a) =
+  a >>> resultById TimingCounterValueId
+
+-- | The counter statistics.
+timingCounterStats :: TimingCounter -> TimingStats
+timingCounterStats (TimingCounter a) =
+  TimingStats (a >>> resultById TimingCounterStatsId)
+
+-- | Represents either finite or infinite queue.
+newtype Queue = Queue ResultTransform
+
+instance ResultTransformer Queue where
+  tr (Queue a) = a
+
+-- | The strategy applied to the enqueueing (input) processes when the finite queue is full.
+enqueueStrategy :: Queue -> ResultTransform
+enqueueStrategy (Queue a) =
+  a >>> resultById EnqueueStrategyId
+
+-- | The strategy applied when storing (in memory) items in the queue.
+enqueueStoringStrategy :: Queue -> ResultTransform
+enqueueStoringStrategy (Queue a) =
+  a >>> resultById EnqueueStoringStrategyId
+
+-- | The strategy applied to the dequeueing (output) processes when the queue is empty.
+dequeueStrategy :: Queue -> ResultTransform
+dequeueStrategy (Queue a) =
+  a >>> resultById DequeueStrategyId
+
+-- | Test whether the queue is empty.
+queueNull :: Queue -> ResultTransform
+queueNull (Queue a) =
+  a >>> resultById QueueNullId
+
+-- | Test whether the finite queue is full.
+queueFull :: Queue -> ResultTransform
+queueFull (Queue a) =
+  a >>> resultById QueueFullId
+
+-- | The finite queue capacity.
+queueMaxCount :: Queue -> ResultTransform
+queueMaxCount (Queue a) =
+  a >>> resultById QueueMaxCountId
+
+-- | Return the current queue size.
+queueCount :: Queue -> ResultTransform
+queueCount (Queue a) =
+  a >>> resultById QueueCountId
+
+-- | Return the queue size statistics.
+queueCountStats :: Queue -> TimingStats
+queueCountStats (Queue a) =
+  TimingStats (a >>> resultById QueueCountStatsId)
+
+-- | Return the total number of input items that were enqueued in the finite queue.
+enqueueCount :: Queue -> ResultTransform
+enqueueCount (Queue a) =
+  a >>> resultById EnqueueCountId
+
+-- | Return the number of lost items for the finite queue.
+enqueueLostCount :: Queue -> ResultTransform
+enqueueLostCount (Queue a) =
+  a >>> resultById EnqueueLostCountId
+
+-- | Return the total number of input items that were stored.
+enqueueStoreCount :: Queue -> ResultTransform
+enqueueStoreCount (Queue a) =
+  a >>> resultById EnqueueStoreCountId
+
+-- | Return the total number of requests for dequeueing the items, not taking
+-- into account the failed attempts to dequeue immediately without suspension.
+dequeueCount :: Queue -> ResultTransform
+dequeueCount (Queue a) =
+  a >>> resultById DequeueCountId
+
+-- | Return the total number of output items that were actually dequeued.
+dequeueExtractCount :: Queue -> ResultTransform
+dequeueExtractCount (Queue a) =
+  a >>> resultById DequeueExtractCountId
+
+-- | Return the load factor: the finite queue size divided by its capacity.
+queueLoadFactor :: Queue -> ResultTransform
+queueLoadFactor (Queue a) =
+  a >>> resultById QueueLoadFactorId
+
+-- | Return the rate of the input items that were enqueued in the finite queue:
+-- how many items per time.
+enqueueRate :: Queue -> ResultTransform
+enqueueRate (Queue a) =
+  a >>> resultById EnqueueRateId
+
+-- | Return the rate of the items that were stored: how many items per time.
+enqueueStoreRate :: Queue -> ResultTransform
+enqueueStoreRate (Queue a) =
+  a >>> resultById EnqueueStoreRateId
+
+-- | Return the rate of the requests for dequeueing the items: how many
+-- requests per time. It does not include the failed attempts to dequeue
+-- immediately without suspension.
+dequeueRate :: Queue -> ResultTransform
+dequeueRate (Queue a) =
+  a >>> resultById DequeueRateId
+
+-- | Return the rate of the output items that were dequeued: how many items per time.
+dequeueExtractRate :: Queue -> ResultTransform
+dequeueExtractRate (Queue a) =
+  a >>> resultById DequeueExtractRateId
+
+-- | Return the wait time from the time at which the item was stored in
+-- the queue to the time at which it was dequeued.
+queueWaitTime :: Queue -> SamplingStats
+queueWaitTime (Queue a) =
+  SamplingStats (a >>> resultById QueueWaitTimeId)
+
+-- | Return the total wait time for the finite queue from the time at which
+-- the enqueueing operation was initiated to the time at which the item was dequeued.
+queueTotalWaitTime :: Queue -> SamplingStats
+queueTotalWaitTime (Queue a) =
+  SamplingStats (a >>> resultById QueueTotalWaitTimeId)
+
+-- | Return the wait time from the time at which the item was stored in
+-- the queue to the time at which it was dequeued.
+enqueueWaitTime :: Queue -> SamplingStats
+enqueueWaitTime (Queue a) =
+  SamplingStats (a >>> resultById EnqueueWaitTimeId)
+
+-- | Return the dequeue wait time from the time at which the item was requested
+-- for dequeueing to the time at which it was actually dequeued.
+dequeueWaitTime :: Queue -> SamplingStats
+dequeueWaitTime (Queue a) =
+  SamplingStats (a >>> resultById DequeueWaitTimeId)
+
+-- | Return a long-term average queue rate calculated as the average queue size
+-- divided by the average wait time.
+queueRate :: Queue -> ResultTransform
+queueRate (Queue a) =
+  a >>> resultById QueueRateId
+
+-- | Accumulates the statistics about that how long the arrived events are processed.
+newtype ArrivalTimer = ArrivalTimer ResultTransform
+
+instance ResultTransformer ArrivalTimer where
+  tr (ArrivalTimer a) = a
+
+-- | Return the statistics about that how long the arrived events were processed.
+arrivalProcessingTime :: ArrivalTimer -> SamplingStats
+arrivalProcessingTime (ArrivalTimer a) =
+  SamplingStats (a >>> resultById ArrivalProcessingTimeId)
+
+-- | It models the server that prodives a service.
+newtype Server = Server ResultTransform
+
+instance ResultTransformer Server where
+  tr (Server a) = a
+
+-- | The initial state of the server.
+serverInitState :: Server -> ResultTransform
+serverInitState (Server a) =
+  a >>> resultById ServerInitStateId
+
+-- | Return the current state of the server.
+serverState :: Server -> ResultTransform
+serverState (Server a) =
+  a >>> resultById ServerStateId
+
+-- | Return the counted total time when the server was locked while
+-- awaiting the input.
+serverTotalInputWaitTime :: Server -> ResultTransform
+serverTotalInputWaitTime (Server a) =
+  a >>> resultById ServerTotalInputWaitTimeId
+
+-- | Return the counted total time spent by the server while
+-- processing the tasks.
+serverTotalProcessingTime :: Server -> ResultTransform
+serverTotalProcessingTime (Server a) =
+  a >>> resultById ServerTotalProcessingTimeId
+
+-- | Return the counted total time when the server was locked while
+-- trying to deliver the output.
+serverTotalOutputWaitTime :: Server -> ResultTransform
+serverTotalOutputWaitTime (Server a) =
+  a >>> resultById ServerTotalOutputWaitTimeId
+
+-- | Return the counted total time spent by the server while it was
+-- preempted waiting for the further proceeding.
+serverTotalPreemptionTime :: Server -> ResultTransform
+serverTotalPreemptionTime (Server a) =
+  a >>> resultById ServerTotalPreemptionTimeId
+
+-- | Return the statistics of the time when the server was locked
+-- while awaiting the input.
+serverInputWaitTime :: Server -> SamplingStats
+serverInputWaitTime (Server a) =
+  SamplingStats (a >>> resultById ServerInputWaitTimeId)
+
+-- | Return the statistics of the time spent by the server while
+-- processing the tasks.
+serverProcessingTime :: Server -> SamplingStats
+serverProcessingTime (Server a) =
+  SamplingStats (a >>> resultById ServerProcessingTimeId)
+
+-- | Return the statistics of the time when the server was locked
+-- while trying to deliver the output.
+serverOutputWaitTime :: Server -> SamplingStats
+serverOutputWaitTime (Server a) =
+  SamplingStats (a >>> resultById ServerOutputWaitTimeId)
+
+-- | Return the statistics of the time spent by the server while
+-- it was preempted waiting for the further proceeding.
+serverPreemptionTime :: Server -> SamplingStats
+serverPreemptionTime (Server a) =
+  SamplingStats (a >>> resultById ServerPreemptionTimeId)
+
+-- | It returns the factor changing from 0 to 1, which estimates
+-- how often the server was awaiting for the next input task.
+serverInputWaitFactor :: Server -> ResultTransform
+serverInputWaitFactor (Server a) =
+  a >>> resultById ServerInputWaitFactorId
+
+-- | It returns the factor changing from 0 to 1, which estimates
+-- how often the server was busy with direct processing its tasks.
+serverProcessingFactor :: Server -> ResultTransform
+serverProcessingFactor (Server a) =
+  a >>> resultById ServerProcessingFactorId
+
+-- | It returns the factor changing from 0 to 1, which estimates
+-- how often the server was locked trying to deliver the output
+-- after the task is finished.
+serverOutputWaitFactor :: Server -> ResultTransform
+serverOutputWaitFactor (Server a) =
+  a >>> resultById ServerOutputWaitFactorId
+
+-- | It returns the factor changing from 0 to 1, which estimates
+-- how often the server was preempted waiting for the further proceeding.
+serverPreemptionFactor :: Server -> ResultTransform
+serverPreemptionFactor (Server a) =
+  a >>> resultById ServerPreemptionFactorId
+
+-- | It models an activity that can be utilised.
+newtype Activity = Activity ResultTransform
+
+instance ResultTransformer Activity where
+  tr (Activity a) = a
+
+-- | The initial state of the activity.
+activityInitState :: Activity -> ResultTransform
+activityInitState (Activity a) =
+  a >>> resultById ActivityInitStateId
+
+-- | Return the current state of the activity.
+activityState :: Activity -> ResultTransform
+activityState (Activity a) =
+  a >>> resultById ActivityStateId
+
+-- | Return the counted total time when the activity was utilised.
+activityTotalUtilisationTime :: Activity -> ResultTransform
+activityTotalUtilisationTime (Activity a) =
+  a >>> resultById ActivityTotalUtilisationTimeId
+
+-- | Return the counted total time when the activity was idle.
+activityTotalIdleTime :: Activity -> ResultTransform
+activityTotalIdleTime (Activity a) =
+  a >>> resultById ActivityTotalIdleTimeId
+
+-- | Return the counted total time when the activity was preemted
+-- waiting for the further proceeding.
+activityTotalPreemptionTime :: Activity -> ResultTransform
+activityTotalPreemptionTime (Activity a) =
+  a >>> resultById ActivityTotalPreemptionTimeId
+
+-- | Return the statistics for the time when the activity was utilised.
+activityUtilisationTime :: Activity -> SamplingStats
+activityUtilisationTime (Activity a) =
+  SamplingStats (a >>> resultById ActivityUtilisationTimeId)
+
+-- | Return the statistics for the time when the activity was idle.
+activityIdleTime :: Activity -> SamplingStats
+activityIdleTime (Activity a) =
+  SamplingStats (a >>> resultById ActivityIdleTimeId)
+
+-- | Return the statistics for the time when the activity was preempted
+-- waiting for the further proceeding.
+activityPreemptionTime :: Activity -> SamplingStats
+activityPreemptionTime (Activity a) =
+  SamplingStats (a >>> resultById ActivityPreemptionTimeId)
+
+-- | It returns the factor changing from 0 to 1, which estimates how often
+-- the activity was utilised.
+activityUtilisationFactor :: Activity -> ResultTransform
+activityUtilisationFactor (Activity a) =
+  a >>> resultById ActivityUtilisationFactorId
+
+-- | It returns the factor changing from 0 to 1, which estimates how often
+-- the activity was idle.
+activityIdleFactor :: Activity -> ResultTransform
+activityIdleFactor (Activity a) =
+  a >>> resultById ActivityIdleFactorId
+
+-- | It returns the factor changing from 0 to 1, which estimates how often
+-- the activity was preempted waiting for the further proceeding.
+activityPreemptionFactor :: Activity -> ResultTransform
+activityPreemptionFactor (Activity a) =
+  a >>> resultById ActivityPreemptionFactorId
+
+-- | The resource which can be acquired and then released.
+newtype Resource = Resource ResultTransform
+
+instance ResultTransformer Resource where
+  tr (Resource a) = a
+
+-- | Return the current available count of the resource.
+resourceCount :: Resource -> ResultTransform
+resourceCount (Resource a) =
+  a >>> resultById ResourceCountId
+
+-- | Return the statistics for the available count of the resource.
+resourceCountStats :: Resource -> TimingStats
+resourceCountStats (Resource a) =
+  TimingStats (a >>> resultById ResourceCountStatsId)
+
+-- | Return the current utilisation count of the resource.
+resourceUtilisationCount :: Resource -> ResultTransform
+resourceUtilisationCount (Resource a) =
+  a >>> resultById ResourceUtilisationCountId
+
+-- | Return the statistics for the utilisation count of the resource.
+resourceUtilisationCountStats :: Resource -> TimingStats
+resourceUtilisationCountStats (Resource a) =
+  TimingStats (a >>> resultById ResourceUtilisationCountStatsId)
+
+-- | Return the current queue length of the resource.
+resourceQueueCount :: Resource -> ResultTransform
+resourceQueueCount (Resource a) =
+  a >>> resultById ResourceQueueCountId
+
+-- | Return the statistics for the queue length of the resource.
+resourceQueueCountStats :: Resource -> TimingStats
+resourceQueueCountStats (Resource a) =
+  TimingStats (a >>> resultById ResourceQueueCountStatsId)
+
+-- | Return the total wait time of the resource.
+resourceTotalWaitTime :: Resource -> ResultTransform
+resourceTotalWaitTime (Resource a) =
+  a >>> resultById ResourceTotalWaitTimeId
+
+-- | Return the statistics for the wait time of the resource.
+resourceWaitTime :: Resource -> SamplingStats
+resourceWaitTime (Resource a) =
+  SamplingStats (a >>> resultById ResourceWaitTimeId)
diff --git a/Simulation/Aivika/Server/Random.hs b/Simulation/Aivika/Server/Random.hs
--- a/Simulation/Aivika/Server/Random.hs
+++ b/Simulation/Aivika/Server/Random.hs
@@ -15,19 +15,32 @@
 module Simulation.Aivika.Server.Random
        (newRandomUniformServer,
         newRandomUniformIntServer,
+        newRandomTriangularServer,
         newRandomNormalServer,
+        newRandomLogNormalServer,
         newRandomExponentialServer,
         newRandomErlangServer,
         newRandomPoissonServer,
         newRandomBinomialServer,
+        newRandomGammaServer,
+        newRandomBetaServer,
+        newRandomWeibullServer,
+        newRandomDiscreteServer,
         newPreemptibleRandomUniformServer,
         newPreemptibleRandomUniformIntServer,
+        newPreemptibleRandomTriangularServer,
         newPreemptibleRandomNormalServer,
+        newPreemptibleRandomLogNormalServer,
         newPreemptibleRandomExponentialServer,
         newPreemptibleRandomErlangServer,
         newPreemptibleRandomPoissonServer,
-        newPreemptibleRandomBinomialServer) where
+        newPreemptibleRandomBinomialServer,
+        newPreemptibleRandomGammaServer,
+        newPreemptibleRandomBetaServer,
+        newPreemptibleRandomWeibullServer,
+        newPreemptibleRandomDiscreteServer) where
 
+import Simulation.Aivika.Generator
 import Simulation.Aivika.Simulation
 import Simulation.Aivika.Process
 import Simulation.Aivika.Process.Random
@@ -62,6 +75,22 @@
   newPreemptibleRandomUniformIntServer False
 
 -- | Create a new server that holds the process for a random time interval
+-- having the triangular distribution, when processing every input element.
+--
+-- By default, it is assumed that the server process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomTriangularServer :: Double
+                             -- ^ the minimum time interval
+                             -> Double
+                             -- ^ the median of the time interval
+                             -> Double
+                             -- ^ the maximum time interval
+                             -> Simulation (Server () a a)
+newRandomTriangularServer =
+  newPreemptibleRandomTriangularServer False
+
+-- | Create a new server that holds the process for a random time interval
 -- distributed normally, when processing every input element.
 --
 -- By default, it is assumed that the server process cannot be preempted,
@@ -76,6 +105,22 @@
   newPreemptibleRandomNormalServer False
          
 -- | Create a new server that holds the process for a random time interval
+-- having the lognormal distribution, when processing every input element.
+--
+-- By default, it is assumed that the server process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomLogNormalServer :: Double
+                            -- ^ the mean of a normal distribution which
+                            -- this distribution is derived from
+                            -> Double
+                            -- ^ the deviation of a normal distribution which
+                            -- this distribution is derived from
+                            -> Simulation (Server () a a)
+newRandomLogNormalServer =
+  newPreemptibleRandomLogNormalServer False
+
+-- | Create a new server that holds the process for a random time interval
 -- distributed exponentially with the specified mean (the reciprocal of the rate),
 -- when processing every input element.
 --
@@ -132,6 +177,63 @@
   newPreemptibleRandomBinomialServer False
 
 -- | Create a new server that holds the process for a random time interval
+-- having the Gamma distribution with the specified shape and scale,
+-- when processing every input element.
+--
+-- By default, it is assumed that the server process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomGammaServer :: Double
+                        -- ^ the shape
+                        -> Double
+                        -- ^ the scale (a reciprocal of the rate)
+                        -> Simulation (Server () a a)
+newRandomGammaServer =
+  newPreemptibleRandomGammaServer False
+
+-- | Create a new server that holds the process for a random time interval
+-- having the Beta distribution with the specified shape parameters (alpha and beta),
+-- when processing every input element.
+--
+-- By default, it is assumed that the server process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomBetaServer :: Double
+                       -- ^ shape (alpha)
+                       -> Double
+                       -- ^ shape (beta)
+                       -> Simulation (Server () a a)
+newRandomBetaServer =
+  newPreemptibleRandomBetaServer False
+
+-- | Create a new server that holds the process for a random time interval
+-- having the Weibull distribution with the specified shape and scale,
+-- when processing every input element.
+--
+-- By default, it is assumed that the server process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomWeibullServer :: Double
+                          -- ^ shape
+                          -> Double
+                          -- ^ scale
+                          -> Simulation (Server () a a)
+newRandomWeibullServer =
+  newPreemptibleRandomWeibullServer False
+
+-- | Create a new server that holds the process for a random time interval
+-- having the specified discrete distribution, when processing every input element.
+--
+-- By default, it is assumed that the server process cannot be preempted,
+-- because the handling of possible task preemption is rather costly
+-- operation.
+newRandomDiscreteServer :: DiscretePDF Double
+                           -- ^ the discrete probability density function
+                           -> Simulation (Server () a a)
+newRandomDiscreteServer =
+  newPreemptibleRandomDiscreteServer False
+
+-- | Create a new server that holds the process for a random time interval
 -- distributed uniformly, when processing every input element.
 newPreemptibleRandomUniformServer :: Bool
                                      -- ^ whether the server process can be preempted
@@ -160,6 +262,22 @@
      return a
 
 -- | Create a new server that holds the process for a random time interval
+-- having the triangular distribution, when processing every input element.
+newPreemptibleRandomTriangularServer :: Bool
+                                        -- ^ whether the server process can be preempted
+                                        -> Double
+                                        -- ^ the minimum time interval
+                                        -> Double
+                                        -- ^ the median of the time interval
+                                        -> Double
+                                        -- ^ the maximum time interval
+                                        -> Simulation (Server () a a)
+newPreemptibleRandomTriangularServer preemptible min median max =
+  newPreemptibleServer preemptible $ \a ->
+  do randomTriangularProcess_ min median max
+     return a
+
+-- | Create a new server that holds the process for a random time interval
 -- distributed normally, when processing every input element.
 newPreemptibleRandomNormalServer :: Bool
                                     -- ^ whether the server process can be preempted
@@ -172,8 +290,24 @@
   newPreemptibleServer preemptible $ \a ->
   do randomNormalProcess_ mu nu
      return a
-         
+
 -- | Create a new server that holds the process for a random time interval
+-- having the lognormal distribution, when processing every input element.
+newPreemptibleRandomLogNormalServer :: Bool
+                                       -- ^ whether the server process can be preempted
+                                       -> Double
+                                       -- ^ the mean of a normal distribution which
+                                       -- this distribution is derived from
+                                       -> Double
+                                       -- ^ the deviation of a normal distribution which
+                                       -- this distribution is derived from
+                                       -> Simulation (Server () a a)
+newPreemptibleRandomLogNormalServer preemptible mu nu =
+  newPreemptibleServer preemptible $ \a ->
+  do randomLogNormalProcess_ mu nu
+     return a
+
+-- | Create a new server that holds the process for a random time interval
 -- distributed exponentially with the specified mean (the reciprocal of the rate),
 -- when processing every input element.
 newPreemptibleRandomExponentialServer :: Bool
@@ -227,4 +361,61 @@
 newPreemptibleRandomBinomialServer preemptible prob trials =
   newPreemptibleServer preemptible $ \a ->
   do randomBinomialProcess_ prob trials
+     return a
+
+-- | Create a new server that holds the process for a random time interval
+-- having the Gamma distribution with the specified shape and scale,
+-- when processing every input element.
+newPreemptibleRandomGammaServer :: Bool
+                                   -- ^ whether the server process can be preempted
+                                   -> Double
+                                   -- ^ the shape
+                                   -> Double
+                                   -- ^ the scale (a reciprocal of the rate)
+                                   -> Simulation (Server () a a)
+newPreemptibleRandomGammaServer preemptible kappa theta =
+  newPreemptibleServer preemptible $ \a ->
+  do randomGammaProcess_ kappa theta
+     return a
+
+-- | Create a new server that holds the process for a random time interval
+-- having the Beta distribution with the specified shape parameters (alpha and beta),
+-- when processing every input element.
+newPreemptibleRandomBetaServer :: Bool
+                                  -- ^ whether the server process can be preempted
+                                  -> Double
+                                  -- ^ shape (alpha)
+                                  -> Double
+                                  -- ^ shape (beta)
+                                  -> Simulation (Server () a a)
+newPreemptibleRandomBetaServer preemptible alpha beta =
+  newPreemptibleServer preemptible $ \a ->
+  do randomBetaProcess_ alpha beta
+     return a
+
+-- | Create a new server that holds the process for a random time interval
+-- having the Weibull distribution with the specified shape and scale,
+-- when processing every input element.
+newPreemptibleRandomWeibullServer :: Bool
+                                  -- ^ whether the server process can be preempted
+                                  -> Double
+                                  -- ^ shape
+                                  -> Double
+                                  -- ^ scale
+                                  -> Simulation (Server () a a)
+newPreemptibleRandomWeibullServer preemptible alpha beta =
+  newPreemptibleServer preemptible $ \a ->
+  do randomWeibullProcess_ alpha beta
+     return a
+
+-- | Create a new server that holds the process for a random time interval
+-- having the specified discrete distribution, when processing every input element.
+newPreemptibleRandomDiscreteServer :: Bool
+                                      -- ^ whether the server process can be preempted
+                                      -> DiscretePDF Double
+                                      -- ^ the discrete probability density function
+                                      -> Simulation (Server () a a)
+newPreemptibleRandomDiscreteServer preemptible dpdf =
+  newPreemptibleServer preemptible $ \a ->
+  do randomDiscreteProcess_ dpdf
      return a
diff --git a/Simulation/Aivika/Stream/Random.hs b/Simulation/Aivika/Stream/Random.hs
--- a/Simulation/Aivika/Stream/Random.hs
+++ b/Simulation/Aivika/Stream/Random.hs
@@ -16,15 +16,22 @@
         randomStream,
         randomUniformStream,
         randomUniformIntStream,
+        randomTriangularStream,
         randomNormalStream,
+        randomLogNormalStream,
         randomExponentialStream,
         randomErlangStream,
         randomPoissonStream,
-        randomBinomialStream) where
+        randomBinomialStream,
+        randomGammaStream,
+        randomBetaStream,
+        randomWeibullStream,
+        randomDiscreteStream) where
 
 import Control.Monad
 import Control.Monad.Trans
 
+import Simulation.Aivika.Generator
 import Simulation.Aivika.Parameter
 import Simulation.Aivika.Parameter.Random
 import Simulation.Aivika.Simulation
@@ -66,7 +73,7 @@
                                    Just t0 -> Just delay }
        return (arrival, Cons $ loop (Just t2))
 
--- | Create a new stream with delays distributed uniformly.
+-- | Create a new stream with random delays distributed uniformly.
 randomUniformStream :: Double
                        -- ^ the minimum delay
                        -> Double
@@ -78,7 +85,7 @@
   randomUniform min max >>= \x ->
   return (x, x)
 
--- | Create a new stream with integer delays distributed uniformly.
+-- | Create a new stream with integer random delays distributed uniformly.
 randomUniformIntStream :: Int
                           -- ^ the minimum delay
                           -> Int
@@ -90,7 +97,21 @@
   randomUniformInt min max >>= \x ->
   return (fromIntegral x, x)
 
--- | Create a new stream with delays distributed normally.
+-- | Create a new stream with random delays having the triangular distribution.
+randomTriangularStream :: Double
+                          -- ^ the minimum delay
+                          -> Double
+                          -- ^ the median of the delay
+                          -> Double
+                          -- ^ the maximum delay
+                          -> Stream (Arrival Double)
+                          -- ^ the stream of random events with the delays generated
+randomTriangularStream min median max =
+  randomStream $
+  randomTriangular min median max >>= \x ->
+  return (x, x)
+
+-- | Create a new stream with random delays distributed normally.
 randomNormalStream :: Double
                       -- ^ the mean delay
                       -> Double
@@ -101,8 +122,22 @@
   randomStream $
   randomNormal mu nu >>= \x ->
   return (x, x)
-         
--- | Return a new stream with delays distibuted exponentially with the specified mean
+
+-- | Create a new stream with random delays having the lognormal distribution.
+randomLogNormalStream :: Double
+                         -- ^ the mean of a normal distribution which
+                         -- this distribution is derived from
+                         -> Double
+                         -- ^ the deviation of a normal distribution which
+                         -- this distribution is derived from
+                         -> Stream (Arrival Double)
+                         -- ^ the stream of random events with the delays generated
+randomLogNormalStream mu nu =
+  randomStream $
+  randomLogNormal mu nu >>= \x ->
+  return (x, x)
+
+-- | Return a new stream with random delays distibuted exponentially with the specified mean
 -- (the reciprocal of the rate).
 randomExponentialStream :: Double
                            -- ^ the mean delay (the reciprocal of the rate)
@@ -113,7 +148,7 @@
   randomExponential mu >>= \x ->
   return (x, x)
          
--- | Return a new stream with delays having the Erlang distribution with the specified
+-- | Return a new stream with random delays having the Erlang distribution with the specified
 -- scale (the reciprocal of the rate) and shape parameters.
 randomErlangStream :: Double
                       -- ^ the scale (the reciprocal of the rate)
@@ -126,7 +161,7 @@
   randomErlang beta m >>= \x ->
   return (x, x)
 
--- | Return a new stream with delays having the Poisson distribution with
+-- | Return a new stream with random delays having the Poisson distribution with
 -- the specified mean.
 randomPoissonStream :: Double
                        -- ^ the mean delay
@@ -137,7 +172,7 @@
   randomPoisson mu >>= \x ->
   return (fromIntegral x, x)
 
--- | Return a new stream with delays having the binomial distribution with the specified
+-- | Return a new stream with random delays having the binomial distribution with the specified
 -- probability and trials.
 randomBinomialStream :: Double
                         -- ^ the probability
@@ -149,3 +184,52 @@
   randomStream $
   randomBinomial prob trials >>= \x ->
   return (fromIntegral x, x)
+
+-- | Return a new stream with random delays having the Gamma distribution by the specified
+-- shape and scale.
+randomGammaStream :: Double
+                     -- ^ the shape
+                     -> Double
+                     -- ^ the scale (a reciprocal of the rate)
+                     -> Stream (Arrival Double)
+                     -- ^ the stream of random events with the delays generated
+randomGammaStream kappa theta =
+  randomStream $
+  randomGamma kappa theta >>= \x ->
+  return (x, x)
+
+-- | Return a new stream with random delays having the Beta distribution by the specified
+-- shape parameters (alpha and beta).
+randomBetaStream :: Double
+                    -- ^ the shape (alpha)
+                    -> Double
+                    -- ^ the shape (beta)
+                    -> Stream (Arrival Double)
+                    -- ^ the stream of random events with the delays generated
+randomBetaStream alpha beta =
+  randomStream $
+  randomBeta alpha beta >>= \x ->
+  return (x, x)
+
+-- | Return a new stream with random delays having the Weibull distribution by the specified
+-- shape and scale.
+randomWeibullStream :: Double
+                       -- ^ shape
+                       -> Double
+                       -- ^ scale
+                       -> Stream (Arrival Double)
+                       -- ^ the stream of random events with the delays generated
+randomWeibullStream alpha beta =
+  randomStream $
+  randomWeibull alpha beta >>= \x ->
+  return (x, x)
+
+-- | Return a new stream with random delays having the specified discrete distribution.
+randomDiscreteStream :: DiscretePDF Double
+                        -- ^ the discrete probability density function
+                        -> Stream (Arrival Double)
+                        -- ^ the stream of random events with the delays generated
+randomDiscreteStream dpdf =
+  randomStream $
+  randomDiscrete dpdf >>= \x ->
+  return (x, x)
diff --git a/Simulation/Aivika/Vector.hs b/Simulation/Aivika/Vector.hs
--- a/Simulation/Aivika/Vector.hs
+++ b/Simulation/Aivika/Vector.hs
@@ -22,7 +22,10 @@
         vectorBinarySearch,
         vectorInsert,
         vectorDeleteAt,
+        vectorDelete,
+        vectorDeleteBy,
         vectorIndex,
+        vectorIndexBy,
         freezeVector) where 
 
 import Data.Array
@@ -183,3 +186,37 @@
                      then return index
                      else loop $ index + 1
      loop 0
+     
+-- | Return an index of the item satisfying the predicate or -1.     
+vectorIndexBy :: Vector a -> (a -> Bool) -> IO Int
+vectorIndexBy vector pred =
+  do count <- readIORef (vectorCountRef vector)
+     array <- readIORef (vectorArrayRef vector)
+     let loop index =
+           if index >= count
+           then return $ -1
+           else do x <- readArray array index
+                   if pred x
+                     then return index
+                     else loop $ index + 1
+     loop 0
+
+-- | Remove the specified element and return a flag indicating
+-- whether the element was found and removed.
+vectorDelete :: Eq a => Vector a -> a -> IO Bool
+vectorDelete vector item =
+  do index <- vectorIndex vector item
+     if index >= 0
+       then do vectorDeleteAt vector index
+               return True
+       else return False
+            
+-- | Remove an element by the specified predicate and return the element if found.
+vectorDeleteBy :: Vector a -> (a -> Bool) -> IO (Maybe a)
+vectorDeleteBy vector pred =
+  do index <- vectorIndexBy vector pred
+     if index >= 0
+       then do a <- readVector vector index
+               vectorDeleteAt vector index
+               return (Just a)
+       else return Nothing
diff --git a/Simulation/Aivika/Vector/Unboxed.hs b/Simulation/Aivika/Vector/Unboxed.hs
--- a/Simulation/Aivika/Vector/Unboxed.hs
+++ b/Simulation/Aivika/Vector/Unboxed.hs
@@ -22,7 +22,10 @@
         vectorBinarySearch,
         vectorInsert,
         vectorDeleteAt,
+        vectorDelete,
+        vectorDeleteBy,
         vectorIndex,
+        vectorIndexBy,
         freezeVector) where 
 
 import Data.Array
@@ -186,3 +189,36 @@
                      else loop $ index + 1
      loop 0
      
+-- | Return an index of the item satisfying the predicate or -1.     
+vectorIndexBy :: Unboxed a => Vector a -> (a -> Bool) -> IO Int
+vectorIndexBy vector pred =
+  do count <- readIORef (vectorCountRef vector)
+     array <- readIORef (vectorArrayRef vector)
+     let loop index =
+           if index >= count
+           then return $ -1
+           else do x <- readArray array index
+                   if pred x
+                     then return index
+                     else loop $ index + 1
+     loop 0
+
+-- | Remove the specified element and return a flag indicating
+-- whether the element was found and removed.
+vectorDelete :: (Unboxed a, Eq a) => Vector a -> a -> IO Bool
+vectorDelete vector item =
+  do index <- vectorIndex vector item
+     if index >= 0
+       then do vectorDeleteAt vector index
+               return True
+       else return False
+            
+-- | Remove an element by the specified predicate and return a flag indicating
+-- whether the element was found and removed.
+vectorDeleteBy :: Unboxed a => Vector a -> (a -> Bool) -> IO Bool
+vectorDeleteBy vector pred =
+  do index <- vectorIndexBy vector pred
+     if index >= 0
+       then do vectorDeleteAt vector index
+               return True
+       else return False
diff --git a/aivika.cabal b/aivika.cabal
--- a/aivika.cabal
+++ b/aivika.cabal
@@ -1,5 +1,5 @@
 name:            aivika
-version:         4.1.1
+version:         4.2
 synopsis:        A multi-paradigm simulation library
 description:
     Aivika is a multi-paradigm simulation library with a strong emphasis
@@ -119,6 +119,8 @@
                      examples/TimeOutWait.hs
                      examples/PingPong.hs
                      examples/PortOperations.hs
+                     examples/SingleLaneTraffic.hs
+                     examples/RenegingFromQueue.hs
                      CHANGELOG.md
 
 library
@@ -137,6 +139,7 @@
                      Simulation.Aivika.Dynamics.Memo.Unboxed
                      Simulation.Aivika.Dynamics.Random
                      Simulation.Aivika.Event
+                     Simulation.Aivika.Gate
                      Simulation.Aivika.Generator
                      Simulation.Aivika.Net
                      Simulation.Aivika.Net.Random
@@ -158,6 +161,7 @@
                      Simulation.Aivika.Resource.Preemption
                      Simulation.Aivika.Resource.Preemption.Base
                      Simulation.Aivika.Results.Locale
+                     Simulation.Aivika.Results.Transform
                      Simulation.Aivika.Results
                      Simulation.Aivika.Results.IO
                      Simulation.Aivika.Server
@@ -209,6 +213,8 @@
                         ExistentialQuantification,
                         TypeFamilies,
                         DeriveDataTypeable,
+                        TypeSynonymInstances,
+                        RankNTypes,
                         CPP
                      
     ghc-options:     -O2
diff --git a/examples/RenegingFromQueue.hs b/examples/RenegingFromQueue.hs
new file mode 100644
--- /dev/null
+++ b/examples/RenegingFromQueue.hs
@@ -0,0 +1,85 @@
+
+-- Example: Reneging from a Queue
+--
+-- It is described in [1]. This is section 7.8.
+-- 
+-- This example models customers arriving to a queue and leaving a queue
+-- after a prescribed period of time. The time between arrivals is 
+-- exponentially distributed with a mean of 20 minutes. The service time 
+-- for customers is uniformly distributed between 15 and 25 minutes. 
+-- Customers will only wait for service if the waiting time is less than
+-- a renege time which is lognormally distributed with a mean of 10 minutes
+-- and a standard deviation of 2 minutes. It is desired to estimate the time
+-- in the system for those customers served, the percent of customers that
+-- renege and the length of the waiting time.
+-- 
+-- [1] A. Alan B. Pritsker, Simulation with Visual SLAM and AweSim, 2nd ed.
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Data.Array
+
+import Simulation.Aivika
+import qualified Simulation.Aivika.Queue.Infinite as IQ
+
+-- | The simulation specs.
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 10000.0,
+                spcDT = 0.1,
+                spcMethod = RungeKutta4,
+                spcGeneratorType = SimpleGenerator }
+
+mu = 10
+nu = 2
+
+logMu = log mu - logNu * logNu / 2
+logNu = sqrt $ log ((nu * nu) / (mu * mu) + 1) 
+
+model :: Simulation Results
+model = do
+  let customers = randomExponentialStream 20
+  timeInSystem <- newRef emptySamplingStats
+  renegeEvents <- newRef (0 :: Int)
+  queue <- runEventInStartTime $ IQ.newFCFSQueue
+  let queueProcess =
+        do x <- IQ.dequeue queue
+           randomUniformProcess_ 15 25
+           t <- liftDynamics time
+           liftEvent $
+             modifyRef timeInSystem $
+             addSamplingStats (t - arrivalTime x)
+           queueProcess
+      renegingProcess x =
+        do randomLogNormalProcess_ logMu logNu
+           liftEvent $
+             do f <- IQ.queueDelete queue x
+                when f $ modifyRef renegeEvents ((+) 1)
+  runProcessInStartTime queueProcess
+  runProcessInStartTime $
+    flip consumeStream customers $ \x ->
+    liftEvent $
+    do IQ.enqueue queue x
+       runProcess $ renegingProcess x
+  return $
+    results
+    [resultSource
+     "queue" "the queue"
+     queue,
+     --
+     resultSource
+     "timeInSystem" "time in system"
+     timeInSystem,
+     --
+     resultSource
+     "renegeEvents" "renege events"
+     renegeEvents]
+
+modelSummary :: Simulation Results
+modelSummary =
+  fmap resultSummary model
+
+main =
+  do printSimulationResultsInStopTime
+       printResultSourceInEnglish
+       modelSummary specs
diff --git a/examples/SingleLaneTraffic.hs b/examples/SingleLaneTraffic.hs
new file mode 100644
--- /dev/null
+++ b/examples/SingleLaneTraffic.hs
@@ -0,0 +1,146 @@
+
+-- Example: Single-Lane Traffic Analysis 
+--
+-- It is described in different sources [1, 2]. So, this is chapter 15 of [2] and section 6.18 of [1].
+--
+-- The system to be modeled in this example consists of the traffic flow from
+-- two directions along a two-lane road, one lane of which has been closed for
+-- 500 meters for repairs. Traffic lights have been placed at each end of
+-- the closed lane to control the flow of traffic through the repair section.
+-- The lights allow traffic to flow for a specified time interval from only
+-- one direction. When a light turns green, the waiting cars start and pass
+-- the light every two seconds. If a car arrives at a green light when there
+-- are no waiting cars, the car passes through the light without delay. The car
+-- arrival pattern is exponentially distributed, with an average of 9 seconds
+-- between cars from direction 1 and 12 seconds between cars from direction 2.
+-- A light cycle consists of green in direction 1, both red, green in direction 2,
+-- both red, and then the cycle is repeated. Both lights remain red for 55 seconds
+-- to allow the cars in transit to leave the repair section before traffic from
+-- the other direction can be initiated.
+-- 
+-- The objective is to simulate the above system to determine values for
+-- the green time for direction 1 and the green time for direction 2 which
+-- yield a low average waiting time for all cars.
+-- 
+-- [1] A. Alan B. Pritsker, Simulation with Visual SLAM and AweSim, 2nd ed.
+-- [2] Труб И.И., Объектно-ориентированное моделирование на C++: Учебный курс. - СПб.: Питер, 2006
+
+import Control.Monad
+import Control.Monad.Trans
+import Control.Arrow
+
+import Data.Array
+
+import Simulation.Aivika
+import qualified Simulation.Aivika.Resource as R
+
+-- | The simulation specs.
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 3600.0,
+                spcDT = 0.1,
+                spcMethod = RungeKutta4,
+                spcGeneratorType = SimpleGenerator }
+
+data LightTime =
+  LightTime { greenLightTime1 :: Double,
+              greenLightTime2 :: Double }
+
+model :: LightTime -> Simulation Results
+model lightTime = do
+  let greenLightTime =
+        array (1, 2)
+        [(1, return $ greenLightTime1 lightTime :: Event Double),
+         (2, return $ greenLightTime2 lightTime :: Event Double)]
+  waitTime1 <- newRef emptySamplingStats
+  waitTime2 <- newRef emptySamplingStats
+  let waitTime =
+        array (1, 2) [(1, waitTime1), (2, waitTime2)]
+  start1 <-
+    runEventInStartTime $
+    R.newFCFSResource 1
+  start2 <-
+    runEventInStartTime $
+    R.newFCFSResource 1
+  let start =
+        array (1, 2) [(1, start1), (2, start2)]
+  light1 <- newGateClosed
+  light2 <- newGateClosed
+  let stream1 = randomExponentialStream 9
+      stream2 = randomExponentialStream 12
+  runProcessInStartTime $
+    flip consumeStream stream1 $ \x ->
+    liftEvent $
+    runProcess $
+    do R.requestResource start1
+       awaitGateOpened light1
+       t <- liftDynamics time
+       liftEvent $
+         modifyRef waitTime1 $
+         addSamplingStats (t - arrivalTime x)
+       when (t > arrivalTime x) $
+         holdProcess 2
+       R.releaseResource start1
+  runProcessInStartTime $
+    flip consumeStream stream2 $ \x ->
+    liftEvent $
+    runProcess $
+    do R.requestResource start2
+       awaitGateOpened light2
+       t <- liftDynamics time
+       liftEvent $
+         modifyRef waitTime2 $
+         addSamplingStats (t - arrivalTime x)
+       when (t > arrivalTime x) $
+         holdProcess 2
+       R.releaseResource start2
+  let lighting =
+        do holdProcess 55
+           liftEvent $
+             openGate light1
+           holdProcess $
+             greenLightTime1 lightTime
+           liftEvent $
+             closeGate light1
+           holdProcess 55
+           liftEvent $
+             openGate light2
+           holdProcess $
+             greenLightTime2 lightTime
+           liftEvent $
+             closeGate light2
+           lighting
+  runProcessInStartTime lighting
+  return $
+    results
+    [resultSource
+     "start" "Start Resource"
+     start,
+     --
+     resultSource
+     "waitTime" "Wait Time"
+     waitTime,
+     --
+     resultSource
+     "greenLightTime" "Green Light Time"
+     greenLightTime]
+
+modelSummary :: LightTime -> Simulation Results
+modelSummary lightTime =
+  fmap resultSummary $ model lightTime
+
+lightTime1 = LightTime 60 45
+lightTime2 = LightTime 80 60
+lightTime3 = LightTime 40 30
+
+model1 = model lightTime1
+model2 = model lightTime2
+model3 = model lightTime3
+
+modelSummary1 = fmap resultSummary model1
+modelSummary2 = fmap resultSummary model2
+modelSummary3 = fmap resultSummary model3
+
+main =
+  do printSimulationResultsInStopTime
+       printResultSourceInEnglish
+       modelSummary1 specs