{-# LANGUAGE TemplateHaskell #-}
module Mcmc.Proposal.Simplex
(
Simplex (toVector),
simplexUniform,
simplexFromVector,
dirichlet,
beta,
)
where
import Data.Aeson
import Data.Aeson.TH
import qualified Data.Vector.Unboxed as V
import Mcmc.Proposal
import Mcmc.Statistics.Types
import Numeric.Log
import Statistics.Distribution
import Statistics.Distribution.Beta
import Statistics.Distribution.Dirichlet
newtype Simplex = SimplexUnsafe {Simplex -> Vector TuningParameter
toVector :: V.Vector Double}
deriving (Simplex -> Simplex -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Simplex -> Simplex -> Bool
$c/= :: Simplex -> Simplex -> Bool
== :: Simplex -> Simplex -> Bool
$c== :: Simplex -> Simplex -> Bool
Eq, Dimension -> Simplex -> ShowS
[Simplex] -> ShowS
Simplex -> String
forall a.
(Dimension -> a -> ShowS)
-> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Simplex] -> ShowS
$cshowList :: [Simplex] -> ShowS
show :: Simplex -> String
$cshow :: Simplex -> String
showsPrec :: Dimension -> Simplex -> ShowS
$cshowsPrec :: Dimension -> Simplex -> ShowS
Show)
$(deriveJSON defaultOptions ''Simplex)
eps :: Double
eps :: TuningParameter
eps = TuningParameter
1e-14
isNormalized :: V.Vector Double -> Bool
isNormalized :: Vector TuningParameter -> Bool
isNormalized Vector TuningParameter
v
| forall a. Num a => a -> a
abs (forall a. (Unbox a, Num a) => Vector a -> a
V.sum Vector TuningParameter
v forall a. Num a => a -> a -> a
- TuningParameter
1.0) forall a. Ord a => a -> a -> Bool
> TuningParameter
eps = Bool
False
| Bool
otherwise = Bool
True
isNegative :: V.Vector Double -> Bool
isNegative :: Vector TuningParameter -> Bool
isNegative = forall a. Unbox a => (a -> Bool) -> Vector a -> Bool
V.any (forall a. Ord a => a -> a -> Bool
< TuningParameter
0)
simplexFromVector :: V.Vector Double -> Either String Simplex
simplexFromVector :: Vector TuningParameter -> Either String Simplex
simplexFromVector Vector TuningParameter
v
| forall a. Unbox a => Vector a -> Bool
V.null Vector TuningParameter
v = forall a b. a -> Either a b
Left String
"simplexFromVector: Vector is empty."
| Vector TuningParameter -> Bool
isNegative Vector TuningParameter
v = forall a b. a -> Either a b
Left String
"simplexFromVector: Vector contains negative elements."
| Bool -> Bool
not (Vector TuningParameter -> Bool
isNormalized Vector TuningParameter
v) = forall a b. a -> Either a b
Left String
"simplexFromVector: Vector is not normalized."
| Bool
otherwise = forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ Vector TuningParameter -> Simplex
SimplexUnsafe Vector TuningParameter
v
simplexUniform :: Dimension -> Simplex
simplexUniform :: Dimension -> Simplex
simplexUniform Dimension
k = forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either forall a. HasCallStack => String -> a
error forall a. a -> a
id forall a b. (a -> b) -> a -> b
$ Vector TuningParameter -> Either String Simplex
simplexFromVector forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Dimension -> a -> Vector a
V.replicate Dimension
k (TuningParameter
1.0 forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Dimension
k)
getTuningFunction :: TuningParameter -> (TuningParameter -> TuningParameter)
getTuningFunction :: TuningParameter -> TuningParameter -> TuningParameter
getTuningFunction TuningParameter
t = (forall a. Fractional a => a -> a -> a
/ TuningParameter
t'')
where
t' :: TuningParameter
t' = TuningParameter
t forall a. Fractional a => a -> a -> a
/ TuningParameter
10000
t'' :: TuningParameter
t'' = forall a. Floating a => a -> a
sqrt TuningParameter
t'
dirichletPFunction :: TuningParameter -> PFunction Simplex
dirichletPFunction :: TuningParameter -> PFunction Simplex
dirichletPFunction TuningParameter
t (SimplexUnsafe Vector TuningParameter
xs) IOGenM StdGen
g = do
let ddXs :: DirichletDistribution
ddXs = forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either forall a. HasCallStack => String -> a
error forall a. a -> a
id forall a b. (a -> b) -> a -> b
$ Vector TuningParameter -> Either String DirichletDistribution
dirichletDistribution forall a b. (a -> b) -> a -> b
$ forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map TuningParameter -> TuningParameter
tf Vector TuningParameter
xs
Vector TuningParameter
ys <- forall g (m :: * -> *).
StatefulGen g m =>
DirichletDistribution -> g -> m (Vector TuningParameter)
dirichletSample DirichletDistribution
ddXs IOGenM StdGen
g
let eitherDdYs :: Either String DirichletDistribution
eitherDdYs = Vector TuningParameter -> Either String DirichletDistribution
dirichletDistribution forall a b. (a -> b) -> a -> b
$ forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map TuningParameter -> TuningParameter
tf Vector TuningParameter
ys
let r :: Jacobian
r = case Either String DirichletDistribution
eitherDdYs of
Left String
_ -> Jacobian
0
Right DirichletDistribution
ddYs -> DirichletDistribution -> Vector TuningParameter -> Jacobian
dirichletDensity DirichletDistribution
ddYs Vector TuningParameter
xs forall a. Fractional a => a -> a -> a
/ DirichletDistribution -> Vector TuningParameter -> Jacobian
dirichletDensity DirichletDistribution
ddXs Vector TuningParameter
ys
forall (f :: * -> *) a. Applicative f => a -> f a
pure (forall a. a -> Jacobian -> Jacobian -> PResult a
Propose (Vector TuningParameter -> Simplex
SimplexUnsafe Vector TuningParameter
ys) Jacobian
r Jacobian
1.0, forall a. Maybe a
Nothing)
where
tf :: TuningParameter -> TuningParameter
tf = TuningParameter -> TuningParameter -> TuningParameter
getTuningFunction TuningParameter
t
dirichlet :: PDimension -> PName -> PWeight -> Tune -> Proposal Simplex
dirichlet :: PDimension -> PName -> PWeight -> Tune -> Proposal Simplex
dirichlet = forall a.
PDescription
-> (TuningParameter -> PFunction a)
-> PSpeed
-> PDimension
-> PName
-> PWeight
-> Tune
-> Proposal a
createProposal (String -> PDescription
PDescription String
"Dirichlet") TuningParameter -> PFunction Simplex
dirichletPFunction PSpeed
PFast
betaPFunction :: Dimension -> TuningParameter -> PFunction Simplex
betaPFunction :: Dimension -> TuningParameter -> PFunction Simplex
betaPFunction Dimension
i TuningParameter
t (SimplexUnsafe Vector TuningParameter
xs) IOGenM StdGen
g = do
let aX :: TuningParameter
aX = TuningParameter
xI
bX :: TuningParameter
bX = TuningParameter
xsSum forall a. Num a => a -> a -> a
- TuningParameter
xI
bdXI :: BetaDistribution
bdXI = TuningParameter -> TuningParameter -> BetaDistribution
betaDistr (TuningParameter -> TuningParameter
tf TuningParameter
aX) (TuningParameter -> TuningParameter
tf TuningParameter
bX)
TuningParameter
yI <- forall d g (m :: * -> *).
(ContGen d, StatefulGen g m) =>
d -> g -> m TuningParameter
genContVar BetaDistribution
bdXI IOGenM StdGen
g
let aY :: TuningParameter
aY = TuningParameter
yI
bY :: TuningParameter
bY = TuningParameter
1.0 forall a. Num a => a -> a -> a
- TuningParameter
yI
eitherBdYI :: Maybe BetaDistribution
eitherBdYI = TuningParameter -> TuningParameter -> Maybe BetaDistribution
betaDistrE (TuningParameter -> TuningParameter
tf TuningParameter
aY) (TuningParameter -> TuningParameter
tf TuningParameter
bY)
let r :: Jacobian
r = case Maybe BetaDistribution
eitherBdYI of
Maybe BetaDistribution
Nothing -> Jacobian
0
Just BetaDistribution
bdYI -> forall a. a -> Log a
Exp forall a b. (a -> b) -> a -> b
$ forall d. ContDistr d => d -> TuningParameter -> TuningParameter
logDensity BetaDistribution
bdYI TuningParameter
xI forall a. Num a => a -> a -> a
- forall d. ContDistr d => d -> TuningParameter -> TuningParameter
logDensity BetaDistribution
bdXI TuningParameter
yI
ja1 :: TuningParameter
ja1 = TuningParameter
bY forall a. Fractional a => a -> a -> a
/ TuningParameter
bX
jac :: Jacobian
jac = forall a. a -> Log a
Exp forall a b. (a -> b) -> a -> b
$ forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall a. Unbox a => Vector a -> Dimension
V.length Vector TuningParameter
xs forall a. Num a => a -> a -> a
- Dimension
2) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log TuningParameter
ja1
let
nf :: TuningParameter -> TuningParameter
nf TuningParameter
x = TuningParameter
x forall a. Num a => a -> a -> a
* TuningParameter
ja1
ys :: Vector TuningParameter
ys = forall a. Unbox a => Dimension -> (Dimension -> a) -> Vector a
V.generate (forall a. Unbox a => Vector a -> Dimension
V.length Vector TuningParameter
xs) (\Dimension
j -> if Dimension
i forall a. Eq a => a -> a -> Bool
== Dimension
j then TuningParameter
yI else TuningParameter -> TuningParameter
nf (Vector TuningParameter
xs forall a. Unbox a => Vector a -> Dimension -> a
V.! Dimension
j))
y :: Simplex
y = forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either forall a. HasCallStack => String -> a
error forall a. a -> a
id forall a b. (a -> b) -> a -> b
$ Vector TuningParameter -> Either String Simplex
simplexFromVector Vector TuningParameter
ys
forall (f :: * -> *) a. Applicative f => a -> f a
pure (forall a. a -> Jacobian -> Jacobian -> PResult a
Propose Simplex
y Jacobian
r Jacobian
jac, forall a. Maybe a
Nothing)
where
xI :: TuningParameter
xI = Vector TuningParameter
xs forall a. Unbox a => Vector a -> Dimension -> a
V.! Dimension
i
xsSum :: TuningParameter
xsSum = forall a. (Unbox a, Num a) => Vector a -> a
V.sum Vector TuningParameter
xs
tf :: TuningParameter -> TuningParameter
tf = TuningParameter -> TuningParameter -> TuningParameter
getTuningFunction TuningParameter
t
beta :: Dimension -> PName -> PWeight -> Tune -> Proposal Simplex
beta :: Dimension -> PName -> PWeight -> Tune -> Proposal Simplex
beta Dimension
i = forall a.
PDescription
-> (TuningParameter -> PFunction a)
-> PSpeed
-> PDimension
-> PName
-> PWeight
-> Tune
-> Proposal a
createProposal PDescription
description (Dimension -> TuningParameter -> PFunction Simplex
betaPFunction Dimension
i) PSpeed
PFast (Dimension -> PDimension
PDimension Dimension
2)
where
description :: PDescription
description = String -> PDescription
PDescription forall a b. (a -> b) -> a -> b
$ String
"Beta; coordinate: " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show Dimension
i