fit-ul-ad.dx

'# Unadjusted Langevin using Dex with AD

-- load some generic utility functions
import djwutils

'## now read and process the data

dat = unsafe_io \. read_file "../pima.data"
AsList(_, tab) = parse_tsv ' ' dat
atab = map (\l. cons "1.0" l) tab
att = map (\r. list2tab r :: (Fin 9)=>String) atab
xStr = map (\r. slice r 0 (Fin 8)) att
xmb = map (\r. map parseString r) xStr :: _=>(Fin 8)=>(Maybe Float)
x = map (\r. map from_just r) xmb :: _=>(Fin 8)=>Float
yStrM = map (\r. slice r 8 (Fin 1)) att
yStr = (transpose yStrM)[0@_]
y = map (\s. select (s == "Yes") 1.0 0.0) yStr

x
y

'## now set up for MCMC

def ll(b: (Fin 8)=>Float) -> Float =
  neg $ sum (log (map (\ x. (exp x) + 1) ((map (\ yi. 1 - 2*yi) y)*(x **. b))))

pscale = [10.0, 1, 1, 1, 1, 1, 1, 1] -- prior SDs
prscale = map (\ x. 1.0/x) pscale

def lprior(b: (Fin 8)=>Float) -> Float =
  bs = b*prscale
  neg $ sum ((log pscale) + (0.5 .* (bs*bs)))

def lpost(b: (Fin 8)=>Float) -> Float =
  (ll b) + (lprior b)

k = new_key 40

def ulKernel(lpi: (Fin n=>Float) -> Float,
    pre: (Fin n)=>Float, dt: Float,
    b: (Fin n)=>Float, k: Key) -> (Fin n)=>Float given (n) =
  -- glpi = grad lpi _
  glpi = \x. grad lpi x
  sdt = sqrt dt
  spre = sqrt pre
  b + ((0.5)*dt) .* (pre*(glpi b)) +
    sdt .* (spre*(randn_vec k))

pre = [100.0,1,1,1,1,1,25,1] -- diagonal pre-conditioner
kern = \b k. ulKernel lpost pre 1.0e-6 b k
init = [-9.0,0,0,0,0,0,0,0]

out = markov_chain init (\s k. step_n 4000 kern s k) 10000 k

mv = meanAndCovariance out
fst mv -- mean
snd mv -- (co)variance matrix

unsafe_io \. write_file "fit-ul-ad.tsv" (to_tsv out)


-- eof