Maybe open an issue to discuss this idea? Might be too far afield for Distributions, but it’s a cool idea that should go somewhere.
— John On Aug 8, 2014, at 4:31 AM, Gabriel Mitchell <[email protected]> wrote: > 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]> 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]> >> 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) >> >
