Skip to content

Instantly share code, notes, and snippets.

@dlakelan
Created January 2, 2023 03:29
Show Gist options
  • Save dlakelan/6204895ecce9b91040745ed5662daff1 to your computer and use it in GitHub Desktop.
Save dlakelan/6204895ecce9b91040745ed5662daff1 to your computer and use it in GitHub Desktop.
Path Integral Approach For Functional Priors
using Pkg
Pkg.activate(".")
Pkg.add(["Turing","Distributions","ApproxFun","QuadGK","StatsPlots","LinearAlgebra","DataFrames"])
using Turing, Distributions, ApproxFun, QuadGK, StatsPlots, LinearAlgebra, DataFrames
# a prior over functions that go up then come down... First let's see how ApproxFun works
Random.set_global_seed!(123456)
f = Fun(Chebyshev(Interval(0..100)),randn(12))
plot(collect(0:100),map(f,0:100))
@model function updown(n)
c ~ MvNormal(repeat([0.0],n),1000.0^2*I(n))
## candidate functions rise from 0 to 1 over 0-20, stay near constant out to 50,
## and decline exponentially ish down to 0.3 out to 100
## this could be some kind of function for say immune response as a function
## of age in years, or whatever
f = Fun(Chebyshev(Interval(0..100)),c)
a = quadgk(x -> ((f(x) - x/20)/.25)^2,0,20) # near linear rise over x 0..20
Turing.@addlogprob!(-a[1])
b = quadgk(x -> ((f(x) - 1.0)/.1)^2,20,50) # near 1 for x in 20-50
Turing.@addlogprob!(-b[1])
c = quadgk(x -> (((f(x) - (exp(-(x-50)/10)+0.2))/.25)^2),50,100) #exponential decay to 0.2
Turing.@addlogprob!(-c[1])
## relatively smooth functions
ddf = Derivative(2)*f
d = quadgk(x -> (ddf(x)^2),0,100)
Turing.@addlogprob!(-d[1]/100.0/.2)
end
s = Turing.sample(updown(10),NUTS(300,.8),300)
df = DataFrame(s)
p = plot(;legend=false)
x = collect(0.0:100.0)
for i in 1:10:300
plot!(x,map(Fun(Chebyshev(Interval(0..100)),Vector(df[i,3:12])),x))
end
display(p)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment