Reputation: 6818
I want to compute K*es
where K
is an Eigen
matrix (dimension pxp
) and es
is a px1
random binary vector with 1s.
For example if p=5
and t=2
a possible es
is [1,0,1,0,0]'
or [0,0,1,1,0]'
and so on...
How do I easily generate es
with Eigen
?
Upvotes: 1
Views: 656
Reputation: 6818
I came up with even a better solution, which is a combination of std::vector
, Egien::Map
and std::shuffle
.
std::vector<int> esv(p,0);
std::fill_n(esv.begin(),t,1);
Eigen::Map<Eigen::VectorXi> es (esv.data(), esv.size());
std::random_device rd;
std::mt19937 g(rd());
std::shuffle(std::begin(esv), std::end(esv), g);
This solution is memory efficient (since Eigen::Map
doesn't copy esv
) and has the big advantage that if we want to permute es
several times (like in this case), then we just need to repeat std::shuffle(std::begin(esv), std::end(esv), g);
Maybe I'm wrong, but this solution seems more elegant and efficient than the previous ones.
Upvotes: 2
Reputation: 9781
When t
is close to p
, Ryan's method need to generate much more than t
random numbers. To avoid this performance degrade, you could solve your original problem
find
t
different numbers from [0, p) that are uniformly distributed
by the following steps
generate t
uniformly distributed random numbers idx[t]
from [0, p-t+1)
sort these numbers idx[t]
idx[i]+i, i=0,...,t-1
are the result
The code:
VectorXi idx(t);
VectorXd es(p);
es.setConstant(0);
for(int i = 0; i < t; ++i) {
idx(i) = int(double(rand()) / RAND_MAX * (p-t+1));
}
std::sort(idx.data(), idx.data() + idx.size());
for(int i = 0; i < t; ++i) {
es(idx(i)+i) = 1.0;
}
Upvotes: 0
Reputation: 313
So you're using Eigen. I'm not sure what matrix type you're using, but I'll go off the class Eigen::MatrixXd
.
What you need to do is:
The following code should do the trick, although you could implement it other ways.
//Your p and t
int p = 5;
int t = 2;
//px1 matrix
MatrixXd es(1, p);
//Initialize the whole 1xp matrix
for (int i = 0; i < p; ++i)
es(1, i) = 0;
//Get a random position in the 1xp matrix from 0-p
for (int i = 0; i < t; ++i)
{
int randPos = rand() % p;
//If the position was already a 1 and not a 0, get a different random position
while (es(1, randPos) == 1)
randPos = rand() % p;
//Change the random position from a 0 to a 1
es(1, randPos) = 1;
}
Upvotes: 0