RandomVariate on multidimensional Bernoulli array

1. "How can I fix this?"

Map RandomVariate at Level {-2}:

SeedRandom[1]
dists2 = Map[BernoulliDistribution, RandomReal[1, {2, 3}], {-1}];

SeedRandom[1];
Map[RandomVariate, dists2, {-2}] 
{{1, 0, 1}, {0, 0, 0}}

You can also use {Depth @ dists2 - 2} as the level specification:

% == Map[RandomVariate, dists2, {Depth @ dists2 - 2}]

True

Why?

Level[dists2, {-1}]
 {0.817389, 0.11142, 0.789526, 0.187803, 0.241361, 0.0657388}
Map[foo, dists2, {-1}]
  {{BernoulliDistribution[foo[0.817389]], 
  BernoulliDistribution[foo[0.11142]], 
  BernoulliDistribution[foo[0.789526]]}, {BernoulliDistribution[
   foo[0.187803]], BernoulliDistribution[foo[0.241361]], 
  BernoulliDistribution[foo[0.0657388]]}}
Level[dists2, {-2}]
 {BernoulliDistribution[0.817389], BernoulliDistribution[0.11142], 
 BernoulliDistribution[0.789526], BernoulliDistribution[0.187803], 
 BernoulliDistribution[0.241361], BernoulliDistribution[0.0657388]}
Map[foo, dists2, {-2}] 
  {{foo[BernoulliDistribution[0.8173894901710712`]], 
  foo[BernoulliDistribution[0.11141961113123644`]], 
  foo[BernoulliDistribution[0.7895259946338515`]]}, {foo[
   BernoulliDistribution[0.18780314670602638`]], 
  foo[BernoulliDistribution[0.24136096745765045`]], 
  foo[BernoulliDistribution[0.06573875950878105`]]}}

2. "why does this not happen in the one-dimensional case?"

If you map RandomVariate at level {-1} the same issue happens in 1D case too:

SeedRandom[1]
dists1 = Map[BernoulliDistribution, RandomReal[1, 3], {-1}];
SeedRandom[1];
Map[RandomVariate, dists1, {-1}] 

enter image description here

The reason your code works for the 1D case without error is that the default level spec for Map is {1} (which is the same as `{Depth@list1 -2}).


An alternative, faster solution to using BernoulliDistribution is the following:

Set up parameters:

dims = {5, 5, 5};

p = RandomReal[1, dims];

Sample:

UnitStep[p - RandomReal[1, dims]]

If you need to use a distribution object for some reason, you could do the following:

Map[RandomVariate@*BernoulliDistribution, p, {-1}]

The reason why your approach did not work was that once you mapped BernoulliDistribution, level {-1} no longer referred to the elements of the 3D array, but to the values in the BernoulliDistribution[...] expressions.