ts-aggregator/project_template/Aic.cpp

#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);
}