Dev (#60)

* update dependencies

* remove junk files

* Reorganizing (Good idea Kusti)

* internals.md

* update README

.gitignore

@@ -5,3 +5,5 @@ deps/deps.jl
 src/sandbox.jl
 Manifest.toml
 Soss.code-workspace
+tracefile*.info
+cache/

.vscode/settings.json

@@ -0,0 +1,3 @@
+{
+    "julia.runLinter": false
+}

Project.toml

@@ -9,16 +9,20 @@ Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DynamicHMC = "bbc10e6e-7c05-544b-b16e-64fede858acb"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+IRTools = "7869d1d1-7146-5819-86e3-90919afe41df"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
+LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogDensityProblems = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
 NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -33,13 +37,16 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 Stheno = "8188c328-b5d6-583d-959b-9690869a5511"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 TransformVariables = "84d833dd-6860-57f9-a1a7-6da5db126cff"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 MonteCarloMeasurements = ">=0.3.5"
 julia = "^1.2"
 
 [extras]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9"
 
 [targets]
-test = ["Test"]
+test = ["Test", "Weave", "BenchmarkTools"]

README.md

@@ -98,7 +98,7 @@ Often these are easier to work with in terms of `particles` (built using [MonteC
 julia> post = dynamicHMC(m(X=truth.X), (y=truth.y,));
 
 julia> particles(post)
-(β = Particles{Float64,1000}[0.561 ± 0.25, 0.755 ± 0.48],)
+(β = Particles{Float64,1000}[0.558 ± 0.25, 0.768 ± 0.49],)
 
 ````
 
@@ -204,11 +204,11 @@ julia> using BenchmarkTools
 
 julia> 
 @btime logpdf($m2(X=X), $truth)
-  809.202 ns (15 allocations: 496 bytes)
+  802.533 ns (16 allocations: 464 bytes)
 -15.84854642585797
 
 julia> @btime logpdf($m2(X=X), $truth, $codegen)
-  321.932 ns (5 allocations: 208 bytes)
+  324.463 ns (5 allocations: 208 bytes)
 -15.848546425857968
 
 ````

docs/internals.md

@@ -0,0 +1,19 @@
+Soss needs the body of a model to be of the form
+```julia
+begin
+    line_1
+    ⋮
+    line_n
+end
+```
+Each line is syntactically translated into a `Statement`. This is an abstract type, with subtypes `Assign` and `Sample`. For example,
+```julia
+x ~ Normal(μ,σ)
+```
+becomes
+```julia
+Sample(:x, :(Normal(μ,σ)))
+```
+Next, all of the `Sample`s are brought together to build a named tuple mapping each `Symbol` to its `Expr`. This becomes the `dists` field for a `Model`.
+
+Because all of this is entirely syntactic, translating into another form only helps when its done on the right side of `~` or `=`. Otherwise we need another way to represent this information.

examples/2019-11-07-demo.jmd

@@ -0,0 +1,226 @@
+Build this file from Julia with
+
+```
+using Weave
+weave("examples/2019-11-07-demo.jmd"
+    , doctype= "github"      # to produce github-friendly markdown
+    , cache=:refresh         # to force the cache to be refreshed
+    , mod=Main               # to evaluate in the Main module
+    , throw_errors=false     # to include any errors in the output (easier to troubleshoot)
+    , args=Dict(:seed => 6)  # to set the random seed
+    )
+```
+
+
+```julia
+using Revise, Soss, Plots, NamedTupleTools
+seed = WEAVE_ARGS[:seed]
+import Random; Random.seed!(seed)
+```
+
+# Make some fake data
+
+```julia
+mt = @model x begin
+    α = 1.0
+    β = 3.0
+    σ = 0.5
+    yhat = α .+ β .* x
+    n = length(x)
+    y ~ For(n) do j
+        Mix([Normal(yhat[j], σ), Normal(yhat[j],8σ)], [0.8,0.2])
+    end
+end;
+
+x = randn(100); 
+truth = rand(mt(x=x));
+```
+
+
+
+
+
+
+```julia
+xx = range(extrema(truth.x)...,length=100)
+scatter(truth.x,truth.y, legend=false, c=1)
+# plot!(xx, truth.α .+ truth.β .* xx, dpi=300,legend=false, lw=3, c=2)
+```
+
+
+```julia 
+m = @model x begin
+    α ~ Normal()
+    β ~ Normal()
+    σ ~ HalfNormal()
+    yhat = α .+ β .* x
+    n = length(x)
+    y ~ For(n) do j
+        Normal(yhat[j], σ)
+    end
+end;
+```
+
+```julia
+m(x=truth.x)
+```
+
+
+```julia
+post = dynamicHMC(m(x=truth.x), (y=truth.y,)) 
+ppost = particles(post)
+```
+
+```julia
+symlogpdf(m) |> foldConstants |> tolatex |> println
+```
+
+```julia
+symlogpdf(m) |> expandSums |> foldConstants |> tolatex |> println
+```
+
+```julia
+using BenchmarkTools
+@btime logpdf($m(x=x),$truth)
+```
+
+```julia
+@btime logpdf($(m(x=x)),$truth, codegen)
+```
+
+
+
+```julia 
+eachplot(xx, ppost.α .+ ppost.β .* xx, lw=3, dpi=300, color=:black)
+scatter!(truth.x,truth.y, legend=false, c=1)
+```
+
+```julia
+
+pred = predictive(m, :α, :β, :σ)
+```
+
+
+```julia
+
+postpred = map(post) do θ 
+    delete(rand(pred(θ)((x=x,))), :n, :x)
+end |> particles
+
+
+pvals = mean.(truth.y .> postpred.y)
+
+
+# PPC vs x
+scatter(truth.x, pvals, legend=false, dpi=300)
+xlabel!("x")
+ylabel!("Bayesian p-value")
+```
+```julia
+
+# 
+# 
+# # PPC vs y
+scatter(truth.y, pvals, legend=false, dpi=300)
+xlabel!("y")
+ylabel!("Bayesian p-value")
+```
+
+
+
+```julia
+
+using AverageShiftedHistograms
+
+o = ash(pvals, rng=0:0.01:1, kernel=Kernels.cosine,m=8)
+plot(o, legend=false,dpi=300)
+xlabel!("Bayesian p-values")
+```
+
+
+```julia
+m2 = @model x begin
+    α ~ Normal()
+    β ~ Normal()
+    σ ~ HalfNormal()
+    yhat = α .+ β .* x
+    νinv ~ HalfNormal()
+    ν = 1/νinv
+    n = length(x)
+    y ~ For(n) do j
+            StudentT(ν,yhat[j],σ)
+        end
+end;
+```
+
+```julia
+post2 = dynamicHMC(m2(x=truth.x), (y=truth.y,))
+ppost2 = particles(post2)
+```
+
+
+```julia
+eachplot(xx, ppost.α .+ ppost.β .* xx, lw=3, dpi=300, color=2)
+eachplot!(xx, ppost2.α .+ ppost2.β .* xx, lw=3, dpi=300, color=:black)
+scatter!(truth.x,truth.y, legend=false, c=1)
+```
+
+```julia
+pred2 = predictive(m2, setdiff(stochastic(m2), [:y])...)
+```
+
+```julia
+post2pred = map(post2) do θ 
+    delete(rand(pred(θ)((x=x,))), :n, :x)
+end |> particles
+
+```
+
+```julia
+pvals2 = mean.(truth.y .> post2pred.y)
+```
+
+
+# PPC vs x
+```julia
+scatter(truth.x, pvals2, legend=false, dpi=300)
+xlabel!("x")
+ylabel!("Bayesian p-value")
+```
+
+
+# PPC vs y
+```julia
+scatter(truth.y, pvals2, legend=false, dpi=300)
+xlabel!("y")
+ylabel!("Bayesian p-value")
+```
+
+```julia
+o = ash(pvals2, rng=0:0.01:1, kernel=Kernels.cosine,m=8)
+plot(o, legend=false,dpi=300)
+xlabel!("Bayesian p-values")
+```
+
+
+
+using Soss
+
+m = @model begin
+    μ ~ Normal() |> iid(2)
+    σ ~ HalfNormal() |> iid(3)
+    x ~ For(1:2,1:3) do i,j
+        Normal(μ[i], σ[j])
+    end
+end;
+
+truth = rand(m())
+
+post = dynamicHMC(m(), (x=truth.x,)) |> particles
+
+pred = predictive(m,:μ,:σ) 
+
+predpost = pred(post) 
+
+rand(predpost)
+```