fit-ul.dx

'# Unadjusted Langevin using Dex

-- 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

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

def glp(b: (Fin 8)=>Float) -> (Fin 8)=>Float =
  glpr = -b*prscale*prscale
  gll = (transpose x) **. (y - (map (\eta. 1.0/(1.0 + eta)) (exp (-x **. b))))
  glpr + gll

k = new_key 42

def ulKernel(glpi: (Fin n=>Float) -> (Fin n)=>Float,
    pre: (Fin n)=>Float, dt: Float,
    b: (Fin n)=>Float, k: Key) -> (Fin n=>Float) given (n) =
  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 glp 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.tsv" (to_tsv out)


-- eof