Reference elements in matrices in Julia while assigning values in other matrices

Viewed 64

I have a vector X whose elements are zeros and ones. I want to create another vector Z of the same size as X where each element of Z is 0 if the corresponding element in X is zero, otherwise it is a random draw from a. uniform distribution. In Matlab I can easily do this by:

n = 1000;
X = randi([0, 1], [1, n]);
Z(X) = rand(); #Here wherever X takes a value of 1, the element of Z is a draw from a uniform distribution. 

I want to implement this in Julia. Is there a cleaner way of doing this instead of using if conditionals. Thanks!!

3 Answers

Here's one way to do it:


julia> n = 1000;

julia> x = rand(Bool, n);

julia> z = zeros(n);

julia> using Distributions

julia> z[x] .= rand.(Uniform(-10, 10));

julia> z
100-element Vector{Float64}:
 -2.6946644136672004
  0.0
  0.0
  ⋮

You can adjust the parameters of the Uniform distribution to what you need, or leave that argument out if the default [0, 1) range is what you need.

The line z[x] .= rand.(Uniform(-10, 10)) uses Julia's logical indexing (same as MATLAB) and broadcasting features - for every x value that is true, the rand call is made and the result assigned to that element of z.

The advantage of using the broadcast (compared to creating rand(Uniform(-10, 10), count(x)) and assigning that to z[x] for eg.) is that the values are directly assigned in-place to their destination in z, and so there's no extra unnecessary memory allocated (as mentioned by @DNF in the comments).

First of all, your Matlab code doesn't work in Matlab, for two reasons: Firstly because logical indices must be boolean, they cannot be 0 and 1. And secondly, because Z(X) = rand() will draw only a single random number and assign it to all the corresponding elements of Z.

Instead, you may want something like this Matlab code:

X = rand(1, n) > 0.5
Z(X) = rand(sum(X), 1)

In Julia you could do

X = rand(Bool, n)
Z = float.(X) # you have to initialize Z
Z[X] .= rand.()

Edit: Here's an alternative with a comprehension, where you don't need to initialize Z:

X = rand(Bool, n)
Z = [x ? float(x) : rand() for x in X]

Technically, what you are sampling from here is a left-censored uniform distribution -- equivalent to the mixture of a Dirac distibution at 0 and Uniform(0, 1). The next release of Distributions.jl will have an implementation for censored, which will remove the need to do any fancy assignment at all:

Z = rand(censored(Uniform(-1.0, 1.0), lower=0.0), N)

where the extent to the left is chosen so that the mixture components have equal weight.

Related