multiplication of ivec with double

[jhjourdan]

The schur multiplication of an ivec with a double is rounded. That is a problem for Bernoulli and Binomial distributions, at least. For example the following program (that simulate a Bernoulli distribution of parameter 0.141), as a completely unwanted behaviour :

#include <iostream>
#include <armadillo>
#include <boost/random.hpp>
#include <cppbugs/cppbugs.hpp>

using namespace arma;
using namespace cppbugs;

int main() {
  ivec x(1);
  x(0) = 0;
  double p = 0.141;

  std::function<void ()> model = [&]() { };

  MCModel<boost::minstd_rand> m(model);
  m.track<Bernoulli>(x).dbern(p);
  m.sample(1e7,1e6,1e4,1);

  int count = 0;
  for(ivec y : m.getNode(x).history)
    count+=y(0);  

  std::cout << "samples: " << m.getNode(x).history.size() << std::endl;
  std::cout << "samples 1: " << count << std::endl;

  return 0;
}

Thanks,

Jacques-Henri Jourdan

[armstrtw]

yes. it's annoying. I think this issue is related to a more general issue in Armadillo in which type promotion is not implemented.

for instance, if you run something like:

#include #include int main() { arma::ivec iv = arma::ivec(2); iv[0] = 3L; iv[1] = 10L; std::cout << log(iv) << std::endl; }

warmstrong@krypton:~/dvl/c++/test/arma.int.log$ ./a.out 1 2

vs R: warmstrong@krypton:~/dvl/c++/test/arma.int.log$ R

log(c(3L,10L)) [1] 1.098612 2.302585

So, arma preserves integers even if you run a calc that returns doubles...

I send this example to Conrad some time ago, but he hasn't addressed it yet as it would require some extensive changes to Armadillo.

I'm sure I could implement my own version of schur that would implement the type promotion, but I haven't gotten around to it yet...

Would you mind sending me your proper email?

And perhaps you can tell me a little about your research too if you have time (and how it relates to cppbugs).

-Whit

On Wed, Jul 11, 2012 at 5:11 AM, jhjourdan [email protected] wrote:

The schur multiplication of an ivec with a double is rounded. That is a problem for Bernoulli and Binomial distributions, at least. For example the following program (that simulate a Bernoulli distribution of parameter 0.141), as a completely unwanted behaviour :

#include <iostream>
#include <armadillo>
#include <boost/random.hpp>
#include <cppbugs/cppbugs.hpp>

using namespace arma;
using namespace cppbugs;

int main() {
  ivec x(1);
  x(0) = 0;
  double p = 0.141;

  std::function<void ()> model = [&]() { };

  MCModel<boost::minstd_rand> m(model);
  m.track<Bernoulli>(x).dbern(p);
  m.sample(1e7,1e6,1e4,1);

  int count = 0;
  for(ivec y : m.getNode(x).history)
    count+=y(0);

  std::cout << "samples: " << m.getNode(x).history.size() << std::endl;
  std::cout << "samples 1: " << count << std::endl;

  return 0;
}

Thanks,

Jacques-Henri Jourdan

---

Reply to this email directly or view it on GitHub: https://github.com/armstrtw/CppBugs/issues/11