Add tests against analytical expressions for first and second moments.

.travis.yml

@@ -2,12 +2,29 @@
 language: julia
 os:
   - linux
-  - osx
+#  - osx
 julia:
   - release
   - nightly
 notifications:
   email: false
+
+# taken from RCall's travis.yml
+before_install:
+  # linux
+  - if [ "$TRAVIS_OS_NAME" == "linux" ]; then sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys E084DAB9; fi
+  - if [ "$TRAVIS_OS_NAME" == "linux" ]; then sudo add-apt-repository -y "deb http://cran.rstudio.com/bin/linux/ubuntu $(lsb_release -s -c)/"; fi
+  - if [ "$TRAVIS_OS_NAME" == "linux" ]; then sudo apt-get update -qq -y; fi
+  - if [ "$TRAVIS_OS_NAME" == "linux" ]; then sudo apt-get install git r-base r-base-dev r-recommended -y; fi
+
+  # osx
+  # faster than using homebrew/science tap
+  # but no permalink to release download
+ # - if [ "$TRAVIS_OS_NAME" == "osx" ]; then wget https://cran.rstudio.com/bin/macosx/ -r --level=1 --accept-regex "R-[0-9.]*.pkg" -nd; fi
+ # - if [ "$TRAVIS_OS_NAME" == "osx" ]; then sudo installer -pkg R-*.pkg -target /; fi
+# # install statmod
+  - wget https://cran.r-project.org/src/contrib/statmod_1.4.29.tar.gz
+  - sudo R CMD INSTALL statmod_1.4.29.tar.gz
 # uncomment the following lines to override the default test script
 #script:
 #  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi

README.md

@@ -21,3 +21,26 @@ to the statmod library.
 
 Note, that `using DistQuads` will start an R process in the background,
 due to a `library(statmod)` call via RCall.
+
+```julia
+julia> using Distributions, DistQuads
+
+julia> bd = Beta(1.4, 5.4)
+Distributions.Beta{Float64}(α=1.4, β=5.4)
+
+julia> dq = DistQuad(bd)
+DistQuads.DistQuad([0.00185197,0.00773162,0.017613,0.0314164,0.0490303,0.0703119,0.095089,0.123161,0.1543,0.188255  …  0.740476,0.779219,0.815718,0.849678,0.880826,0.90891,0.933707,0.955022,0.972705,0.986694],[0.00387997,0.013685,0.0273551,0.0429761,0.0586422,0.0726105,0.0834732,0.0902888,0.0926473,0.0906605  …  0.00117028,0.000555584,0.000239261,9.17104e-5,3.04716e-5,8.44413e-6,1.83878e-6,2.84652e-7,2.58301e-8,8.47448e-10],Distributions.Beta{Float64}(α=1.4, β=5.4))
+
+julia> mean(dq)
+0.20588235294117635
+
+julia> mean(bd)
+0.20588235294117643
+
+julia> var(dq)
+0.020960873036997594
+
+julia> var(bd)
+0.020960873036997597
+
+```

REQUIRE

@@ -1,3 +1,4 @@
 julia 0.5
 Distributions
 RCall
+Reexport

deps/build.jl

@@ -0,0 +1,17 @@
+using RCall
+try
+    R"options(repos=structure(c(CRAN=\"https://cloud.r-project.org/\")))"
+    R"""pkgTest <- function(x)
+      {
+        if (!require(x,character.only = TRUE, quiet = TRUE))
+        {
+          install.packages(x,dep=TRUE)
+            if(!require(x,character.only = TRUE)) stop(\"Package not found\")
+        }
+      }"""
+    R"pkgTest(\"statmod\")"
+    R"did_it_work <- \"statmod\" %in% loadedNamespaces()"
+    @assert @rget did_it_work
+catch m
+    error("statmod was not installed, build.jl has failed to set up DistQuads.jl. Please report this to the issue tracker https://github.com/pkofod/DistQuads.jl/issues")
+end

src/DistQuads.jl

@@ -1,45 +1,52 @@
 module DistQuads # DistQuads
 
-using Distributions,
-      RCall
-export DistQuad, E
+using RCall,
+      Reexport
+@reexport using Distributions
+export DistQuad, E, mean, var
+import Base: mean, var
 # FIXME Should install statmod via build!
 R"library(statmod)"
 
 type DistQuad
     x
     w
+    d
 end
 nodes(dq::DistQuad) = dq.x
 weights(dq::DistQuad) = dq.w
+distribution(dq::DistQuad) = dq.d
 
-DistQuad(d) = DistQuad(d, 32)
+DistQuad(d; N = 32) = DistQuad(d, N)
 function DistQuad(d::Distributions.Beta, N)
     R"nodesweight<-gauss.quad.prob($N, dist=\"beta\", alpha=$(d.α), beta=$(d.β))"
     nw = @rget nodesweight
-    DistQuad(nw[:nodes], nw[:weights])
+    DistQuad(nw[:nodes], nw[:weights], d)
 end
 
 function DistQuad(d::Distributions.Gamma, N)
-    R"nodesweight<-gauss.quad.prob($N, dist=\"gamma\", alpha=$(d.α), beta=$(d.Θ))"
+    R"nodesweight<-gauss.quad.prob($N, dist=\"gamma\", alpha=$(d.α), beta=$(d.θ))"
     nw = @rget nodesweight
-    DistQuad(nw[:nodes], nw[:weights])
+    DistQuad(nw[:nodes], nw[:weights], d)
 end
 
 function DistQuad(d::Distributions.Normal, N)
     R"nodesweight<-gauss.quad.prob($N, dist=\"normal\", mu=$(d.μ), sigma=$(d.σ))"
     nw = @rget nodesweight
-    DistQuad(nw[:nodes], nw[:weights])
+    DistQuad(nw[:nodes], nw[:weights], d)
 end
 
 function DistQuad(d::Distributions.LogNormal, N)
     R"nodesweight<-gauss.quad.prob($N, dist=\"normal\", mu=$(d.μ), sigma=$(d.σ))"
     nw = @rget nodesweight
-    DistQuad(exp.(nw[:nodes]), nw[:weights])
+    DistQuad(exp.(nw[:nodes]), nw[:weights], d)
 end
 
 
 E(f, dq::DistQuad) = dot(f.(nodes(dq)), weights(dq))
 E(f, d) = E(f, DistQuad(d))
 
+mean(dq::DistQuad) = E(identity, dq)
+var(dq::DistQuad) = E(x->(x-mean(dq))^2, dq)
+
 end # module

test/runtests.jl

@@ -1,5 +1,62 @@
-using DistQuads
+using DistQuads, Distributions
 using Base.Test
 
-# write your own tests here
-@test 1 == 2
+
+@testset "analytic vs quadrature" begin
+    @testset "Beta" begin
+        for a = 0.1:0.1:3.0, b = 0.1:0.1:3.0
+            bd = Beta(a, b)
+            dq = DistQuad(bd)
+            # Test that mean function didn't accidentally regress
+            @test mean(dq) == E(identity, dq)
+            # Test nodes and weights against Distributions
+            @test norm(mean(dq)-mean(bd), Inf) < 1e-8
+            @test norm(var(dq)-var(bd), Inf) < 1e-8
+            # Test nodes and weights against know mean (Distributions could in principle regress)
+            @test norm(mean(dq)-a/(a+b), Inf) < 1e-8
+            @test norm(var(dq)-a*b/((a+b+1)*(a+b)^2), Inf) < 1e-8
+        end
+    end
+    @testset "Gamma" begin
+        for a = 0.1:0.1:3.0, b = 0.1:0.1:3.0
+            bd = Gamma(a, b)
+            dq = DistQuad(bd)
+            # Test that mean function didn't accidentally regress
+            @test mean(dq) == E(identity, dq)
+            # Test nodes and weights against Distributions
+            @test norm(mean(dq)-mean(bd), Inf) < 1e-8
+            @test norm(var(dq)-var(bd), Inf) < 1e-8
+            # Test nodes and weights against know mean (Distributions could in principle regress)
+            @test norm(mean(dq)-a*b, Inf) < 1e-8
+            @test norm(var(dq)-a*b^2, Inf) < 1e-8
+        end
+    end
+    @testset "Normal" begin
+        for μ = 0.1:0.1:3.0, σ = 0.1:0.1:3.0
+            bd = Normal(μ, σ)
+            dq = DistQuad(bd)
+            # Test that mean function didn't accidentally regress
+            @test mean(dq) == E(identity, dq)
+            # Test nodes and weights against Distributions
+            @test norm(mean(dq)-mean(bd), Inf) < 1e-8
+            @test norm(var(dq)-var(bd), Inf) < 1e-8
+            # Test nodes and weights against know mean (Distributions could in principle regress)
+            @test norm(mean(dq)-μ, Inf) < 1e-8
+            @test norm(var(dq)-σ^2, Inf) < 1e-8
+        end
+    end
+    @testset "LogNormal" begin
+        for μ = 0.1:0.1:3.0, σ = 0.1:0.1:2.0
+            bd = LogNormal(μ, σ)
+            dq = DistQuad(bd)
+            # Test that mean function didn't accidentally regress
+            @test mean(dq) == E(identity, dq)
+            # Test nodes and weights against Distributions
+            @test norm(mean(dq)-mean(bd), Inf) < 1e-8
+            @test norm(var(dq)-var(bd), Inf) < 1e-8 # not too precise!
+            # Test nodes and weights against know mean (Distributions could in principle regress)
+            @test norm(mean(dq)-exp(μ+σ^2/2), Inf) < 1e-8
+            @test norm(var(dq)-(exp(σ^2)-1)*exp(2μ+σ^2), Inf) < 1e-8 # not too precise!
+        end
+    end
+end