hello world example of sampling

Project.toml

@@ -16,6 +16,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 SymbolicNumericIntegration = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

examples/sample/dep_bits.jl

@@ -0,0 +1,37 @@
+using Revise
+using Dice
+# using Random
+using Distributions
+
+NSAMPLES=1000
+samples = rand(Bernoulli(0.5), NSAMPLES)
+samples = map(x -> if x == 1 true else false end, samples)
+
+pws = Dict{}(0 => 0.4, 1 => 0.6)
+
+function dice_program(s)
+    param = if s 2/3 else 3/4 end
+    # param = 0.7
+    code = @dice begin
+            flip1 = flip(param)
+            # flip2 = if flip1 flip(4/7) else flip(2/3) end
+            flip2 = s
+            w = DistUInt{2}([flip1, flip2])
+            observe(w < DistUInt{2}([flip(0.5), flip(0.5)]))
+            w
+    end
+    inside_expectation = expectation(code)
+    prob_evidence = pr(allobservations(code))
+    # @show prob_evidence
+    return inside_expectation, prob_evidence[true], pws[s]
+end
+
+
+mean = map(x -> dice_program(x), samples)
+
+num = sum(map(x -> x[1]*x[2]*x[3]/0.5, mean))
+den = sum(map(x -> x[2]*x[3]/0.5, mean))
+print(num/den)
+
+# true distribution ---- 0 -> 0.1, 1 -> 0.2, 2 -> 0.3, 3 -> 0.4
+# true expectation ---- 1

examples/sample/gt.jl

@@ -0,0 +1,30 @@
+using Dice
+# using Random
+using Distributions
+
+NSAMPLES=100000
+samples = rand(Bernoulli(0.5), NSAMPLES)
+samples = map(x -> if x == 1 true else false end, samples)
+
+function dice_program(s)
+    code = @dice begin
+            w = DistUInt{2}([flip(0.5), s])
+            observe(w < DistUInt{2}([flip(0.5), flip(0.5)]))
+            w
+    end
+    inside_expectation = expectation(code)
+    prob_evidence = pr(allobservations(code))
+    # @show prob_evidence
+    return inside_expectation, prob_evidence[true]
+end
+
+
+mean = map(x -> dice_program(x), samples)
+
+num = sum(map(x -> x[1]*x[2], mean))
+den = sum(map(x -> x[2], mean))
+print(num/den)
+
+# true distribution ---- 0 -> 0.5, 1 -> 0.33, 2 -> 0.166, 3 -> 0
+# true expectation ---- 0.66
+

examples/sample/normal_lower_bits_fixed.jl

@@ -0,0 +1,19 @@
+using Revise
+using Dice
+using Distributions
+using Plots
+
+a = bitblast(DistFix{9, 0}, Normal(0, 1), 16, -8.0, 8.0)
+code = @dice begin
+            # observe(prob_equals(a.mantissa.number.bits[5:7], [false, false, false]))
+            b = bitblast(DistFix{9, 0}, Normal(0, 1), 16, -8.0, 8.0)
+            observe(prob_equals(a + b, DistFix{9, 0}(1.0)))
+            a
+end
+
+# p = pr(a.mantissa.number.bits[5:7])
+
+# plot(pr(a))
+
+p = expectation(code)
+plot(p)