I agree that if the intended use is limited to getting a random value one
time then Johan's solution looks really nice. However, if you wanted to
make the same call several times throughout your code and/or to compute
other statistical quantities of interest it might be worth thinking about
something along the lines of
type GenericIndependentProduct <: DiscreteMultivariateDistribution
dim::Integer
v::Vector
GenericIndependentProduct(dim,v) = dim != length(v) ? error("dim does not
match length(s)") : new(dim,v)
end
GenericIndependentProduct(v::Vector) = GenericIndependentProduct(length(v),v
)
rand(X::GenericIndependentProduct) = [rand(X_i) for X_i in X.v]
rand(X::GenericIndependentProduct,N::Integer) = [rand(X) for i in 1:N]
#mean(X::GenericIndependentProduct) = [mean(X_i) for X_i in X.v]
#and so on...DRY suggests a bit of metaprogramming could help
This assumes that it makes sense to think of the samples as part of a
single (now multivariate). One can then write things like
X = GenericIndependentProduct([Bernoulli(p) for p in rand(7)])
x = rand(X,2)
This might seem overengineered, but (at slight risk of getting off topic) I
am compelled to point out that the admittedly terribly named
GenericIndependentProduct is supposed to suggest a relationship to the the
already established scheme for multivariate normals where we have
GenericMvNormal{PDMat},GenericMvNormal{PDiagMat} and
GenericMvNormal{ScalMat} corresponding to arbitrarily correlated,
independent but not necessarily identical and iid normal distributions,
respectively. Both these multivariate normal distributions the distribution
from which Sam L wants to sample are special cases of the probability
models assumed by linear independent component analysis
(http://en.wikipedia.org/wiki/Independent_component_analysis), inasmuch as
there exists a deterministic (linear) transformation of the random
variables that results in a vector of independent scalar random variables.
Is there any sense for if/when/how the Distributions package would like to
represent factorization structure and/or deterministic transformations of
random variables?
On Tuesday, January 14, 2014 12:00:25 AM UTC+1, John Myles White wrote:
>
> FWIW, I really like Johan's approach.
>
> For me, the main difficulty with emulating R's interface is that the
> current version of `rand` is a method that takes exactly one distribution
> object and draws samples from that object. R doesn't have to think about
> the number of distributions objects being passed in, so it exposes a very
> different interface.
>
> One thing we could do is allow users to pass in a vector of distribution
> objects, each of which would contribute one more samples. But then you'd
> end up writing something like rand([Bernoulli(p) for p in [0.1, 0.2,
> 0.3]]). I think that's a little strange as an idiom, but it would mix
> naturally with our existing tools.
>
> -- John
>
> On Jan 13, 2014, at 1:05 PM, Kevin Squire <[email protected]
> <javascript:>> wrote:
>
> Hi Sam, I would suggest submitting an issue against the Distributions.jl
> package.
>
> Cheers, Kevin
>
>
>
>
> On Mon, Jan 13, 2014 at 12:00 PM, Johan Sigfrids <[email protected]
> <javascript:>> wrote:
>
>> I would think using a list comprehension would be a easy and Julian way
>> of doing it:
>>
>> using Distributions
>> p = [0.1, 0.2, 0.3, 0.4, 0.5]
>> [rand(Bernoulli(x)) for x in p]
>>
>>
>>
>> On Monday, January 13, 2014 9:44:02 PM UTC+2, Sam L wrote:
>>>
>>> In my case I'm trying to generate an array of independent Bernoulli
>>> r.v.s where the parameter for each Bernoulli is given in an array p. In R,
>>> I'd do rbinom(length(p), 1, p). It looks like the Distributions
>>> package does not support this sort of thing. I have a few solutions to my
>>> particular problem shown at the bottom of this message but I'm curious
>>> about the best way to do this in general.
>>>
>>> Is this something the Distributions package will support in the future
>>> or is there a better/built in way to do this? It seems to me that it would
>>> come up frequently and having to write a function to vectorize it every
>>> time is a pain.
>>>
>>> Sam
>>>
>>> #Some solutions
>>> rbern(p) = rand(Bernoulli(p))
>>> @vectorize_1arg Real rbern
>>> type RandBernoulli <: Functor{1} end
>>> NumericExtensions.evaluate(::RandBernoulli, p) = rbern(p) #then map
>>>
>>> rbern2(p) = int(rand(length(p)) .< p)
>>>
>>
>
>
