https://github.com/cran/robCompositions
Tip revision: 6cf109eab116e889a3e3bcc1309cbdcc254895e8 authored by Matthias Templ on 25 August 2023, 15:30:06 UTC
version 2.4.1
version 2.4.1
Tip revision: 6cf109e
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;
}