#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
double dBETA(double dY, double dAlpha, double dBeta, bool bLog=false) {
double dLPDF = (dAlpha - 1.0)*log(dY) + (dBeta - 1.0)*log(1.0 - dY) + Rf_lgammafn(dAlpha + dBeta) -
Rf_lgammafn(dAlpha) - Rf_lgammafn(dBeta);
if(!bLog) dLPDF=exp(dLPDF);
return dLPDF;
}
double pBETA(double dY, double dAlpha, double dBeta) {
double dP = Rf_pbeta(dY, dAlpha, dBeta, 1, 0);
return dP;
}
double qBETA(double dP, double dAlpha, double dBeta) {
double dQ = Rf_qbeta(dP,dAlpha,dBeta,1,0);
return dQ;
}
//
double rBETA(double dAlpha, double dBeta){
double dY = Rf_rbeta(dAlpha, dBeta);
return dY;
}
//
arma::vec mBETA(double dAlpha, double dBeta){
arma::vec vMoments(4);
vMoments(0) = dAlpha/(dAlpha + dBeta);
vMoments(1) = dAlpha*dBeta/( pow(dAlpha + dBeta, 2.0)*(dAlpha + dBeta + 1.0));
vMoments(2) = (2.0*(dBeta - dAlpha)*pow(dAlpha + dBeta + 1.0, 0.5))/((dAlpha + dBeta + 2.0)*pow(dAlpha*dBeta,0.5));
vMoments(3) = 6.0*( pow(dAlpha - dBeta, 2.0) * (dAlpha + dBeta + 1.0) - dAlpha*dBeta*(dAlpha + dBeta + 2.0) )/(dAlpha*dBeta*(dAlpha + dBeta + 2.0)*(dAlpha + dBeta + 3.0)) + 3.0;
return vMoments;
}
//
arma::vec beta_Score(double dY, arma::vec vTheta){
double dAlpha = vTheta(0);
double dBeta = vTheta(1);
arma::vec vScore(2);
double dAlpha_s = log(dY) + Rf_digamma(dAlpha + dBeta) - Rf_digamma(dAlpha);
double dBeta_s = log(1.0 - dY) + Rf_digamma(dAlpha + dBeta) - Rf_digamma(dBeta);
vScore(0) = dAlpha_s;
vScore(1) = dBeta_s;
return vScore;
}
arma::mat beta_IM( arma::vec vTheta){
double dAlpha = vTheta(0);
double dBeta = vTheta(1);
arma::mat mIM = zeros(2,2);
double dTrig = Rf_trigamma(dAlpha + dBeta);
mIM(0,0) = Rf_trigamma(dAlpha) - dTrig;
mIM(1,1) = Rf_trigamma(dBeta) - dTrig;
mIM(1,0) = - dTrig;
return mIM;
}
//
//