Revision a53065a09c3fce65a63e137deb5bccb6162e6cff authored by Matthias Templ on 18 November 2020, 20:10 UTC, committed by cran-robot on 18 November 2020, 20:10 UTC
1 parent 9ae1e67
zeros.cpp
``````#include <math.h>
#include <Eigen/Dense>
#include <iostream>
#include <iterator>
#include <vector>
#include <algorithm>
#include <functional>
#include <numeric>
#include "zeros.h"

std::vector<double>
help::divide
(const std::vector<double> & vect, const double & D)
{
assert(D>0 && " Error: dividing by zero..");
std::vector<double> out;
for(auto it:vect) out.push_back(it/D);
return out;
};

double
help::sum
(const std::vector<double> & vect)
{
return std::accumulate(vect.begin(),vect.end(),0.0);
};

std::vector<double>
help::uniform
(const unsigned int & n)
{
std::vector<double> tmp;
tmp.insert(tmp.begin(),n,1.0/(double)(n));
return tmp;
};

double
help::geom_mean
(const std::vector<double> & vect)
{
double out = 1.0;
for(const auto it:vect){
out=out*it;
}
return pow(out, 1.0/vect.size());
};

void
BM
(std::vector<double> & numbers,
const Eigen::Block<Eigen::Map<Eigen::Matrix<double, -1, -1>, 0, Eigen::Stride<0, 0> >, 1, -1, false> & data,
PRIOR p)
{
if(p == PRIOR::DEFAULT)
{
if( data.sum() > (double) (data.size()*data.size()) )  // n > D^2
p = PRIOR::SQ;
else
p = PRIOR::BAYES_LAPLACE;
}

double s;

switch(p){
case PRIOR::PERKS :
{
s = 1;
BM(numbers, data, s);
break;
}
case PRIOR::JEFFREYS :
{
s = (double)(data.size())/2;
BM(numbers, data, s);
break;
}
case PRIOR::BAYES_LAPLACE :
{
s = (double)(data.size());
BM(numbers, data, s);
break;
}
case PRIOR::SQ :
{
s = sqrt(data.sum());
BM(numbers, data, s);
break;
}
default: {}
}

return;
};

void
BM
(std::vector<double> & numbers,
const Eigen::Block<Eigen::Map<Eigen::Matrix<double, -1, -1>, 0, Eigen::Stride<0, 0> >, 1, -1, false> & data,
const double & s, const bool is_strength_inverse)
{
std::vector<double> t = help::uniform(data.size());

assert(s>=0 && " Error (BM): strength must be >=0..");
assert(data.size()==t.size() && " Error(BM): different sizes of input..");

double n = data.sum();
double tmp = 0.0;

double t_tot = 0.0;

for(std::size_t i = 0; i < data.size(); i++)
{   //computing term of the summation
assert(data[i]>=0 && " Error (BM): input must be >=0..");
assert(t[i]>=0 && " Error (BM): input must be >=0..");

if(data(i)==0)
tmp += t[i];

t_tot += t[i];
}

// assert()  NOTE: need to check that t_tot is "equal" to 1
for(std::size_t i = 0; i < data.size(); i++)
{  //applying BM method
if(data[i] == 0)
{
if(is_strength_inverse == false)
numbers.push_back(t[i]*s/(n + s));
else
numbers.push_back(t[i]/(s*n + 1));
}
else
{
if(is_strength_inverse == false)
numbers.push_back(data[i]*(1 - s*tmp/(n+s))/n);
else
numbers.push_back(data[i]*(1 - tmp/(s*n+1))/n);
}
}
return;
}
``````

Computing file changes ...