/// @file power_regression_approximation.cpp
/// @author Vo Hoang Anh, <spyofgame200@gmail.com>
/// @brief This is a simple implementation of the power regression approximation
/// With given array (x[]) and (y[]) at same size (n)
/// It will find constant (a), (b) that provides
/// Approximation function f(x) = a * x[i] ^ b ~ y[i]
///
/// @license Public domain.
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
/// Size limitation
const int LIM = 1e6 + 16;
/// Given array
int n;
long double x[LIM], y[LIM];
/// Constant need to find
long double a, b;
/// set it to (-1) if you only need the approximation function
/// Number of digits after decimal points of (a) and (b)
int w = 15;
/// First min(k, n) number to be comparing f(x[i]) vs (y[i])
int k = 0;
/// Funtion: a * x^b ~ y
long double f(long double x) { return a * pow(x, b); }
int main()
{
/// Input array
cin >> n;
for (int i = 1; i <= n; ++i) cin >> x[i];
for (int i = 1; i <= n; ++i) cin >> y[i];
/// Detail information
cin >> w >> k;
/// Precalculating
long double sigma_ln_x = 0;
long double sigma_ln_y = 0;
long double sigma_ln_xx = 0;
long double sigma_ln_xy = 0;
for (int i = 1; i <= n; ++i)
{
long double tx = log(x[i]);
long double ty = log(y[i]);
sigma_ln_x += tx;
sigma_ln_y += ty;
sigma_ln_xx += tx * tx;
sigma_ln_xy += tx * ty;
}
/// Power approximation constant
b = ((n * sigma_ln_xy) - (sigma_ln_x * sigma_ln_y))
/ ((n * sigma_ln_xx) - (sigma_ln_x * sigma_ln_x));
/// Base approximation constant
a = exp((sigma_ln_y - b * sigma_ln_x) / n);
/// Output such function
cout << setprecision(w) << "Approximation Function: f(x) = " << a << " * x ^(" << b << ")" << endl;
/// Incase you dont need error
if (w == -1) return 0;
/// Calculating error
long double Sigma_error_value_down = 0;
long double Sigma_error_value_up = 0;
long double Sigma_error_percent_down = 0;
long double Sigma_error_percent_up = 0;
long double Sigma_error_value = 0;
long double Sigma_error_percent = 0;
k = min(k, n);
if (k) cout << "\nComparing f(x[i]) vs y[i]:\n";
for (int i = 1; i <= n; ++i)
{
long double Error_value = y[i] - f(x[i]);
long double Error_percent = Error_value / y[i] * 100;
if (i <= k) /// Comparing f(x[i]) vs (y[i])
{
cout << "At: " << setw(6) << i << " | ";
cout << setprecision(w) << fixed << "<" << f(x[i]) << " -> " << y[i] << "> = ";
cout << setprecision(w) << fixed << Error_value << " ~ " << Error_percent << " %" << endl;
}
if (Error_value < 0) Sigma_error_value_down += Error_value;
if (Error_value > 0) Sigma_error_value_up += Error_value;
Sigma_error_value += fabs(Error_value);
if (Error_percent < 0) Sigma_error_percent_down += Error_percent;
if (Error_percent > 0) Sigma_error_percent_up += Error_percent;
Sigma_error_percent += fabs(Error_percent);
}
/// Average error
long double Average_error_value_down = Sigma_error_value_down / n;
long double Average_error_value_up = Sigma_error_value_up / n;
long double Average_error_value = Sigma_error_value / n;
long double Average_error_percent_down = Sigma_error_percent_down / n;
long double Average_error_percent_up = Sigma_error_percent_up / n;
long double Average_error_percent = Sigma_error_percent / n;
/// Output error
cout << "\nFunction Errors: " << setprecision(4) << Average_error_percent << "%\n";
cout << "Sigma_error_value_down " << " = " << setprecision(w) << fixed << Sigma_error_value_down << endl;
cout << "Sigma_error_value_up " << " = " << setprecision(w) << fixed << Sigma_error_value_up << endl;
cout << "Sigma_error_value " << " = " << setprecision(w) << fixed << Sigma_error_value << endl;
cout << "Average_error_value_down " << " = " << setprecision(w) << fixed << Average_error_value_down << endl;
cout << "Average_error_value_up " << " = " << setprecision(w) << fixed << Average_error_value_up << endl;
cout << "Average_error_value " << " = " << setprecision(w) << fixed << Average_error_value << endl;
cout << "Sigma_error_percent_down " << " = " << setprecision(w) << fixed << Sigma_error_percent_down << "%" << endl;
cout << "Sigma_error_percent_up " << " = " << setprecision(w) << fixed << Sigma_error_percent_up << "%" << endl;
cout << "Sigma_error_percent " << " = " << setprecision(w) << fixed << Sigma_error_percent << "%" << endl;
cout << "Average_error_percent_down " << " = " << setprecision(w) << fixed << Average_error_percent_down << "%" << endl;
cout << "Average_error_percent_up " << " = " << setprecision(w) << fixed << Average_error_percent_up << "%" << endl;
cout << "Average_error_percent " << " = " << setprecision(w) << fixed << Average_error_percent << "%" << endl;
return 0;
}