# Create a custom Random.Sampler or Matrix sampler via Distribution.Sampleable

**URL:** https://discourse.julialang.org/t/create-a-custom-random-sampler-or-matrix-sampler-via-distribution-sampleable/79834
**Category:** General Usage
**Tags:** random
**Created:** [April 22, 2022, 12:30pm UTC](https://discourse.julialang.org/t/create-a-custom-random-sampler-or-matrix-sampler-via-distribution-sampleable/79834 "2022-04-22T12:30:14Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![dmetivie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmetivie/32/6926_2.png) [@dmetivie](https://discourse.julialang.org/u/dmetivie)
#### Post date: [April 22, 2022, 12:30pm UTC](https://discourse.julialang.org/t/create-a-custom-random-sampler-or-matrix-sampler-via-distribution-sampleable/79834/1 "2022-04-22T12:30:14Z")

</div>

I am trying to define `rand` methods to “randomize” points `x` represented in a matrix `(n, s)` (`n` is the number of point and `s` is the dimension).  
I could directly define

```julia
function rand(RM::RandomizationMethod)
     x = RM.x
     return shift(x)
end

```

where `shift(x)` can be a random translation e.g. `shift(x) = x+rand(size(x))` and `RM` is a structure containing information on the original set of point to randomize (here directly the points but it can be their bit representation in some base and other properties).  
All my different randomization methods are gathered in a type `RandomizationMethod`.

My questions:

- I have [read](https://docs.julialang.org/en/v1/stdlib/Random/#Hooking-into-the-Random-API) about `Random.Sampler` but I am not sure to understand what they actually do (especially for nontrivial cases) and why one need them.  
In my case are they necessary?  
From what I understood, if I manage to define the appropriate `Sampler` then `rand` and all associate functions `rand!` or `rand(RM, 5)` (for 5 randomization so a `(n,s,5)` array) or `rand(RM, 5,2)` (for a `(n,s,5,2)` array) should work almost directly.  
In particular, I could have a `Sampler` for one randomization and a slightly different if I want to generate multiple ones more efficiently.

- I saw [Matrix Sampler](https://juliastats.org/Distributions.jl/stable/extends/#Matrix-variate-Sampler) from the `Distribution.jl` package. It seems to be doing sort of what I want with the `Sampleable` type.  
Is it more appropriate for my problem?

- I am a bit lost on how to properly use that for matrix sampling: especially the part to specialize my code for multiple realization.  
Here is a sort of MWE of my issue

```julia
using Random
using Distributions: Matrixvariate, Continuous
import Distributions: Sampleable
import Random: _rand!, rand!

abstract type SamplerShifting{F<:Matrixvariate,S<:Continuous} <: Sampleable{F,S} end
Base.size(sampler::SamplerShifting) = size(sampler.points) # the size of each matrix sample
struct Shifter{F<:AbstractMatrix{<:Real}} <: SamplerShifting{Matrixvariate,Continuous}
    points::F
    #... potentially other things
end

function shift!(rng, sampler, x)
    for i in eachindex(x)
        x[i] = sampler.points[i] + rand(rng)
    end
end

function _rand!(rng::AbstractRNG, sampler::SamplerShifting, A::DenseMatrix{T}) where {T<:Real}
    # ... generate a single matrix sample to A
    #... potentially some initialization
    shift!(rng, sampler, A)
    return A #? is it not bad to return the value here for performance ? 
end

function rand!(rng::AbstractRNG, s::Sampleable{Matrixvariate}, A::AbstractMatrix)
    size(A) == size(s) ||
        throw(DimensionMismatch("Output size inconsistent with sample length."))
    _rand!(rng, s, A)
end

m = 5
n = 2^m
s = 10
x = zeros(n, s)
sp = Shifter(x)

rand(sp, 5, 2) #! this work! Returns a Matrix of my matrix samples

#* Now I want to generate mutliple samples to do all the initialization only once
#* so I try the following (which do not work)
function _rand!(rng::AbstractRNG, sampler::SamplerShifting, A::Array{DenseMatrix{T}}) where {T<:Real}
    # ... generate a multiple matrix sample to A
    #... potentially some initialization
    for i in eachindex(A)
        shift!(rng, sampler, view(A, i)) #? this in place replacment do not work proprely
    end
    return A
end

function rand!(rng::AbstractRNG, s::Sampleable{Matrixvariate}, A::Array{DenseMatrix{T}}) where {T<:Real}
    size(A[firstindex(A)]) == size(s) ||
        throw(DimensionMismatch("Output size inconsistent with sample length."))
    _rand!(rng, s, A)
end

rand(sp, 5, 2) #!! Run but this does NOT use the new method for multiple call 
A = [x, x]
rand!(sp, A) #! Run but produce the same realization

```