diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index ecd3e8a..27a99a3 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -1,7 +1,17 @@
 name: CI
 on:
-  - push
-  - pull_request
+  push:
+    branches:
+      - master
+    paths-ignore:
+      - 'LICENSE'
+      - 'README.md'
+      - '.github/workflows/TagBot.yml'
+  pull_request:
+    paths-ignore:
+      - 'LICENSE'
+      - 'README.md'
+      - '.github/workflows/TagBot.yml'
 jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
diff --git a/Project.toml b/Project.toml
index 256b8b1..e3efb5a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.9.3"
+version = "0.9.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -22,15 +22,20 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[weakdeps]
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+
+[extensions]
+MultivariateOrthogonalPolynomialsStatsBaseExt = "StatsBase"
 
 [compat]
 ArrayLayouts = "1.11"
 BandedMatrices = "1"
 BlockArrays = "1.0"
 BlockBandedMatrices = "0.13"
-ClassicalOrthogonalPolynomials = "0.15"
-ContinuumArrays = "0.19.5"
+ClassicalOrthogonalPolynomials = "0.15.8"
+ContinuumArrays = "0.20"
 DomainSets = "0.7"
 FastTransforms = "0.17"
 FillArrays = "1.0"
@@ -39,17 +44,21 @@ InfiniteArrays = "0.15"
 InfiniteLinearAlgebra = "0.9, 0.10"
 LazyArrays = "2.3.1"
 LazyBandedMatrices = "0.11.3"
-QuasiArrays = "0.12.2"
+QuasiArrays = "0.13"
+Random = "1"
 RecurrenceRelationships = "0.2"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"
+StatsBase = "0.34"
 julia = "1.10"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "ForwardDiff", "RecipesBase", "Test"]
+test = ["Base64", "ForwardDiff", "Random", "RecipesBase", "StatsBase", "Test"]
diff --git a/ext/MultivariateOrthogonalPolynomialsStatsBaseExt.jl b/ext/MultivariateOrthogonalPolynomialsStatsBaseExt.jl
new file mode 100644
index 0000000..abb76c9
--- /dev/null
+++ b/ext/MultivariateOrthogonalPolynomialsStatsBaseExt.jl
@@ -0,0 +1,18 @@
+module MultivariateOrthogonalPolynomialsStatsBaseExt
+using MultivariateOrthogonalPolynomials, StatsBase
+using MultivariateOrthogonalPolynomials.StaticArrays
+using MultivariateOrthogonalPolynomials.QuasiArrays
+using MultivariateOrthogonalPolynomials: ExpansionLayout, KronOPLayout
+import MultivariateOrthogonalPolynomials.QuasiArrays: sample_layout
+
+function sample_layout(::ExpansionLayout{KronOPLayout{2}}, f::AbstractQuasiVector, n::Integer)
+    F = reshape(f)
+    x = sample(sum(F; dims=2), n)
+    # x = sample(F[:,y]) # TODO: this should work
+    y = [sample(F[x,:]) for x in x]
+    map(SVector, x, y)
+end
+
+sample_layout(lay::ExpansionLayout{KronOPLayout{2}}, f::AbstractQuasiVector) = only(sample_layout(lay, f, 1))
+
+end
\ No newline at end of file
diff --git a/src/rect.jl b/src/rect.jl
index f1c56fc..cee393f 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -202,4 +202,4 @@ end
 function vec_layout(::KronExpansionLayout{OPLayout, OPLayout}, f)
     A,C,Bc = arguments(f)
     RectPolynomial(A, Bc') * DiagTrav(C)
-end
\ No newline at end of file
+end
diff --git a/test/test_rect.jl b/test/test_rect.jl
index 1136e78..1dd3d0b 100644
--- a/test/test_rect.jl
+++ b/test/test_rect.jl
@@ -1,9 +1,11 @@
-using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, BlockArrays, FillArrays, Base64, LazyBandedMatrices, ArrayLayouts, Test
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, BlockArrays, FillArrays, Base64, LazyBandedMatrices, ArrayLayouts, Random, StatsBase, Test
 using ClassicalOrthogonalPolynomials: expand, coefficients, recurrencecoefficients
 using MultivariateOrthogonalPolynomials: weaklaplacian, ClenshawKron
 using ContinuumArrays: plotgridvalues, ExpansionLayout
 using Base: oneto
 
+Random.seed!(3242)
+
 @testset "RectPolynomial" begin
     @testset "Evaluation" begin
         T = ChebyshevT()
@@ -161,7 +163,7 @@ using Base: oneto
         f = expand(P, splat((x,y) -> 1))
         @test diff(f,(1,0))[SVector(0.1,0.2)] == diff(f,(0,1))[SVector(0.1,0.2)] == 0.0
         f = expand(P, splat((x,y) -> x))
-        @test diff(f,(1,0))[SVector(0.1,0.2)] == 1.0
+        @test diff(f,(1,0))[SVector(0.1,0.2)] ≈ 1.0
         @test diff(f,(0,1))[SVector(0.1,0.2)] == 0.0
         f = expand(P, splat((x,y) -> cos(x*exp(y))))
         @test diff(f,(1,0))[SVector(0.1,0.2)] ≈ -sin(0.1*exp(0.2))*exp(0.2)
@@ -247,5 +249,9 @@ using Base: oneto
         F = reshape(f)
         @test sum(F; dims=1)[1,0.2] ≈ 2.346737615950585
         @test sum(F; dims=2)[0.1,1] ≈ 2.1748993079723618
+        
+        x,y = coordinates(P)
+        @test sample(f) isa SVector
+        @test sum(sample(f, 100_000))/100_000 ≈ [sum(x .* f)/sum(f),sum(y .* f)/sum(f)] rtol=1E-1
     end
 end