{-# LANGUAGE RankNTypes #-}

-- |
-- Module      :  Mcmc.Proposal.Scale
-- Description :  Multiplicative proposals
-- Copyright   :  2021 Dominik Schrempf
-- License     :  GPL-3.0-or-later
--
-- Maintainer  :  dominik.schrempf@gmail.com
-- Stability   :  unstable
-- Portability :  portable
--
-- Creation date: Thu May 14 21:49:23 2020.
module Mcmc.Proposal.Scale
  ( scale,
    scaleUnbiased,
    scaleContrarily,
  )
where

import Mcmc.Proposal
import Mcmc.Proposal.Generic
import Mcmc.Statistics.Types
import Numeric.Log
import Statistics.Distribution.Gamma

-- The actual proposal with tuning parameter. The tuning parameter does not
-- change the mean.
scalePFunction :: Shape Double -> Scale Double -> TuningParameter -> PFunction Double
scalePFunction :: Shape Double
-> Shape Double -> Shape Double -> PFunction (Shape Double)
scalePFunction Shape Double
k Shape Double
th Shape Double
t =
  GammaDistribution
-> (Shape Double -> Shape Double -> Shape Double)
-> Maybe (Shape Double -> Shape Double)
-> Maybe (Shape Double -> Shape Double -> Jacobian)
-> PFunction (Shape Double)
forall d a.
(ContDistr d, ContGen d) =>
d
-> (a -> Shape Double -> a)
-> Maybe (Shape Double -> Shape Double)
-> Maybe (a -> Shape Double -> Jacobian)
-> PFunction a
genericContinuous
    (Shape Double -> Shape Double -> GammaDistribution
gammaDistr (Shape Double
k Shape Double -> Shape Double -> Shape Double
forall a. Fractional a => a -> a -> a
/ Shape Double
t) (Shape Double
th Shape Double -> Shape Double -> Shape Double
forall a. Num a => a -> a -> a
* Shape Double
t))
    Shape Double -> Shape Double -> Shape Double
forall a. Num a => a -> a -> a
(*)
    ((Shape Double -> Shape Double)
-> Maybe (Shape Double -> Shape Double)
forall a. a -> Maybe a
Just Shape Double -> Shape Double
forall a. Fractional a => a -> a
recip)
    ((Shape Double -> Shape Double -> Jacobian)
-> Maybe (Shape Double -> Shape Double -> Jacobian)
forall a. a -> Maybe a
Just Shape Double -> Shape Double -> Jacobian
forall a p. Floating a => p -> a -> Log a
jac)
  where
    jac :: p -> a -> Log a
jac p
_ = a -> Log a
forall a. a -> Log a
Exp (a -> Log a) -> (a -> a) -> a -> Log a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Floating a => a -> a
log (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Fractional a => a -> a
recip

-- | Multiplicative proposal.
--
-- The gamma distribution is used to sample the multiplier. Therefore, this and
-- all derived proposals are log-additive in that they do not change the sign of
-- the state. Further, the value zero is never proposed when having a strictly
-- positive value.
--
-- Consider using 'Mcmc.Proposal.Slide.slide' to allow proposition of values
-- having opposite sign.
scale ::
  Shape Double ->
  Scale Double ->
  PName ->
  PWeight ->
  Tune ->
  Proposal Double
scale :: Shape Double
-> Shape Double
-> PName
-> PWeight
-> Tune
-> Proposal (Shape Double)
scale Shape Double
k Shape Double
th = PDescription
-> (Shape Double -> PFunction (Shape Double))
-> PSpeed
-> PDimension
-> PName
-> PWeight
-> Tune
-> Proposal (Shape Double)
forall a.
PDescription
-> (Shape Double -> PFunction a)
-> PSpeed
-> PDimension
-> PName
-> PWeight
-> Tune
-> Proposal a
createProposal PDescription
description (Shape Double
-> Shape Double -> Shape Double -> PFunction (Shape Double)
scalePFunction Shape Double
k Shape Double
th) PSpeed
PFast (Int -> PDimension
PDimension Int
1)
  where
    description :: PDescription
description = String -> PDescription
PDescription (String -> PDescription) -> String -> PDescription
forall a b. (a -> b) -> a -> b
$ String
"Scale; shape: " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Shape Double -> String
forall a. Show a => a -> String
show Shape Double
k String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
", scale: " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Shape Double -> String
forall a. Show a => a -> String
show Shape Double
th

-- | See 'scale'.
--
-- The scale of the gamma distribution is set to (shape)^{-1}, so that the mean
-- of the gamma distribution is 1.0.
scaleUnbiased ::
  Shape Double ->
  PName ->
  PWeight ->
  Tune ->
  Proposal Double
scaleUnbiased :: Shape Double -> PName -> PWeight -> Tune -> Proposal (Shape Double)
scaleUnbiased Shape Double
k = PDescription
-> (Shape Double -> PFunction (Shape Double))
-> PSpeed
-> PDimension
-> PName
-> PWeight
-> Tune
-> Proposal (Shape Double)
forall a.
PDescription
-> (Shape Double -> PFunction a)
-> PSpeed
-> PDimension
-> PName
-> PWeight
-> Tune
-> Proposal a
createProposal PDescription
description (Shape Double
-> Shape Double -> Shape Double -> PFunction (Shape Double)
scalePFunction Shape Double
k (Shape Double
1 Shape Double -> Shape Double -> Shape Double
forall a. Fractional a => a -> a -> a
/ Shape Double
k)) PSpeed
PFast (Int -> PDimension
PDimension Int
1)
  where
    description :: PDescription
description = String -> PDescription
PDescription (String -> PDescription) -> String -> PDescription
forall a b. (a -> b) -> a -> b
$ String
"Scale unbiased; shape: " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Shape Double -> String
forall a. Show a => a -> String
show Shape Double
k

scaleContrarilyPFunction ::
  Shape Double ->
  Scale Double ->
  TuningParameter ->
  PFunction (Double, Double)
scaleContrarilyPFunction :: Shape Double
-> Shape Double
-> Shape Double
-> PFunction (Shape Double, Shape Double)
scaleContrarilyPFunction Shape Double
k Shape Double
th Shape Double
t =
  GammaDistribution
-> ((Shape Double, Shape Double)
    -> Shape Double -> (Shape Double, Shape Double))
-> Maybe (Shape Double -> Shape Double)
-> Maybe ((Shape Double, Shape Double) -> Shape Double -> Jacobian)
-> PFunction (Shape Double, Shape Double)
forall d a.
(ContDistr d, ContGen d) =>
d
-> (a -> Shape Double -> a)
-> Maybe (Shape Double -> Shape Double)
-> Maybe (a -> Shape Double -> Jacobian)
-> PFunction a
genericContinuous
    (Shape Double -> Shape Double -> GammaDistribution
gammaDistr (Shape Double
k Shape Double -> Shape Double -> Shape Double
forall a. Fractional a => a -> a -> a
/ Shape Double
t) (Shape Double
th Shape Double -> Shape Double -> Shape Double
forall a. Num a => a -> a -> a
* Shape Double
t))
    (Shape Double, Shape Double)
-> Shape Double -> (Shape Double, Shape Double)
forall b. Fractional b => (b, b) -> b -> (b, b)
contra
    ((Shape Double -> Shape Double)
-> Maybe (Shape Double -> Shape Double)
forall a. a -> Maybe a
Just Shape Double -> Shape Double
forall a. Fractional a => a -> a
recip)
    (((Shape Double, Shape Double) -> Shape Double -> Jacobian)
-> Maybe ((Shape Double, Shape Double) -> Shape Double -> Jacobian)
forall a. a -> Maybe a
Just (Shape Double, Shape Double) -> Shape Double -> Jacobian
forall a p. Floating a => p -> a -> Log a
jac)
  where
    contra :: (b, b) -> b -> (b, b)
contra (b
x, b
y) b
u = (b
x b -> b -> b
forall a. Num a => a -> a -> a
* b
u, b
y b -> b -> b
forall a. Fractional a => a -> a -> a
/ b
u)
    jac :: p -> a -> Log a
jac p
_ a
u = a -> Log a
forall a. a -> Log a
Exp (a -> Log a) -> a -> Log a
forall a b. (a -> b) -> a -> b
$ a -> a
forall a. Floating a => a -> a
log (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
$ a -> a
forall a. Fractional a => a -> a
recip (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
$ a
u a -> a -> a
forall a. Num a => a -> a -> a
* a
u

-- | See 'scale'.
--
-- The two values are scaled contrarily so that their product stays constant.
-- Contrary proposals are useful when parameters are confounded.
scaleContrarily ::
  Shape Double ->
  Scale Double ->
  PName ->
  PWeight ->
  Tune ->
  Proposal (Double, Double)
scaleContrarily :: Shape Double
-> Shape Double
-> PName
-> PWeight
-> Tune
-> Proposal (Shape Double, Shape Double)
scaleContrarily Shape Double
k Shape Double
th = PDescription
-> (Shape Double -> PFunction (Shape Double, Shape Double))
-> PSpeed
-> PDimension
-> PName
-> PWeight
-> Tune
-> Proposal (Shape Double, Shape Double)
forall a.
PDescription
-> (Shape Double -> PFunction a)
-> PSpeed
-> PDimension
-> PName
-> PWeight
-> Tune
-> Proposal a
createProposal PDescription
description (Shape Double
-> Shape Double
-> Shape Double
-> PFunction (Shape Double, Shape Double)
scaleContrarilyPFunction Shape Double
k Shape Double
th) 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
"Scale contrariliy; shape: " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Shape Double -> String
forall a. Show a => a -> String
show Shape Double
k String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
", scale: " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Shape Double -> String
forall a. Show a => a -> String
show Shape Double
th