ts-aggregator/project_template/Aic.cpp
2022-12-13 12:36:06 +04:00

203 lines
5.8 KiB
C++
Raw Permalink Blame History

This file contains invisible Unicode characters

This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#include "StdAfx.h"
#include <iostream>
#include <math.h>
#include "Aic.h"
# define M_PI 3.14159265358979323846 /* pi */
using namespace std;
// критерий Акаике. Наследуется от класса Estimation.
// реализует метод "получить значение критерия"
Aic::Aic() {
}
Aic::~Aic(){
}
/*
paramsCount - count of model params,
количество параетров модели
tsLen - lenght of forecasted estimated time seria,
длина оцениваемого спрогнозированного временного ряда
isSecondOrder - AIC second order modification,
модифицированный АИК второго порядка
original - values of the predicted time series,
известные значения прогнозируемого временного ряда
model - forecasted values of the predicted time series,
прогнозированные значения прогнозируемого временного ряда
*/
//Надо перенести в интерфейс модели, логично если модели будут в последствии
//разных типов
//Функция правдоподобия для не линейной модели
double Aic::logLikNLS(vector<double> forecast, vector<double> weights){
int n = forecast.size();
if (weights.size() == 0){
for (int i = 0; i < n; i++){
weights.push_back(1);
}
}
//summ log(w)
double wls = 0;
//wetghted forecasts
double wfor = 0;
for (int i = 0; i < n; i++){
//zero weights hack
if (weights[i] != 0){
wls += log(weights[i]);
}
else{
wls += log((double)1);
}
wfor += weights[i] * forecast[i] * forecast[i];
}
double res = -n * (log((double)2 * M_PI) + 1 - log((double)n) - wls + log((double)wfor)) / 2;
return res;
}
//Функция правдоподобия для линейной модели
double Aic::logLikLM(vector<double> forecast, vector<double> weights){
int n = forecast.size();
if (weights.size() == 0){
for (int i = 0; i < n; i++){
weights.push_back(1);
}
}
else{
int i = 0;
while (i < n){
if (weights[i] == 0){
forecast.erase(forecast.begin() + i);
weights.erase(weights.begin() + i);
n = forecast.size();
}
i++;
}
}
//summ log(w)
double wls = 0;
//wetghted forecasts
double wfor = 0;
for (int i = 0; i < n; i++){
wls += log(weights[i]);
wfor += weights[i] * forecast[i] * forecast[i];
}
double res = 0.5 * (wls - n * (log((double)2 * M_PI) + 1 - log((double)n) + log((double)wfor)));
return res;
}
//Свойство класса модели
int Aic::getParamsCount(vector<double> model){
#pragma message ("Need solve where is it must do")
cout << "log-likelihood function calculate";
return 0;
}
/*
models with special methods:
glm, multinom
+- maxlikeFit, merMod, unmarkedFit, vglm
*/
double Aic::getBicValue(int paramsCount,
vector<double> original, vector<double> model){
int n = (original.size() > model.size() ? model.size() : original.size());
double res = 0;
double lik = 0;
double lik2 = 0;
vector<double> e;
cout << "Bicc " << endl;
cout << " e sum(lik) sum(lik2)" << endl;
for (int i = 0; i < n; i++){
//adderror
e.push_back(original[i] - model[i]);
cout << " " << e[i];
lik += e[i] * e[i];
cout << " " << lik;
lik2 += log(fabs(model[i]));
cout << " " << lik2 << endl;
// MULT lik += (original[i] - model[i])/ model[i];
}
lik = n * log(lik);
//MULT lik +=2 *lik2;
res = lik + log((double)n) * paramsCount;
cout << " " << res << " " << endl;
return res;
}
double Aic::getAiccValue(int paramsCount,
vector<double> original, vector<double> model){
int n = (original.size() > model.size() ? model.size() : original.size());
double res = 0;
double lik = 0;
double lik2 = 0;
vector<double> e;
cout << "Aicc " << endl;
cout << " e sum(lik) sum(lik2)" << endl;
for (int i = 0; i < n; i++){
//adderror
e.push_back(original[i] - model[i]);
cout << " " << e[i];
lik += e[i] * e[i];
cout << " " << lik;
lik2 += log(fabs(model[i]));
cout << " " << lik2 << endl;
// MULT lik += (original[i] - model[i])/ model[i];
}
lik = n * log(lik);
//MULT lik +=2 *lik2;
res = lik + (2 * paramsCount* n)/( n - paramsCount - 1);
cout << " " << res << " " << endl;
return res;
}
/*
valide for next models by package:
aov, clm, clmm, coxme, coxph, gls, lm, lme, lmeik, mer, merMod,
multinom, nlme, nls, polr, rlm, zeroinfl
*/
double Aic::getValue(int paramsCount, // bool isSecondOrder,
vector<double> original, vector<double> model) {
//Validate(original, model);
int n = (original.size() > model.size() ? model.size() : original.size());
////// R default package
//double logLik = Aic::logLikNLS(model, RSSweights(original, model));//Aic::logLikLM(model, RSSweights(original, model));
////R-package version
//if (true){//isSecondOrder){
// return -2 * logLik + 2 * paramsCount;
//}
//else
//{
// return -2 * logLik + 2 * paramsCount *(n / (n - paramsCount - 1));
//}
//lm version from paper
/*double rss = RSS(original, model);
double res = 2 * paramsCount + n * log(rss / (n - 2));
return res;*/
///forecast 5.6 R package
double res = 0;
double lik = 0;
double lik2 = 0;
vector<double> e;
//cout << "Aic " << endl;
//cout << " e sum(lik) sum(lik2)" << endl;
for (int i = 0; i < n; i++){
//adderror
e.push_back(original[i] - model[i]);
//cout << " " << e[i];
lik += e[i] * e[i];
//cout << " " << lik;
lik2 += log(fabs(model[i]));
//cout << " " << lik2 << endl;
// MULT lik += (original[i] - model[i])/ model[i];
}
lik = n * log(lik);
//MULT lik +=2 *lik2;
res = lik + 2 * paramsCount;
//cout << " " << res << " " << endl;
return res;
}
double Aic::getValue(vector<double> original, vector<double> model) {
return this->getValue(4, original, model);
}