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