Skip to content

Instantly share code, notes, and snippets.

@newpolaris

newpolaris/datagen_plot.py

Last active May 3, 2021 03:08
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save newpolaris/5da6f82cad0e6c9767fac83969162913 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save newpolaris/5da6f82cad0e6c9767fac83969162913 to your computer and use it in GitHub Desktop.
Download ZIP
polynomialRegiression
Raw
datagen_plot.py
// https://himbopsa.tistory.com/13
import numpy as np
import matplotlib.pyplot as plt
x = 5 * np.random.rand(1,150)
y = -2 * pow(x,3) + 9 *pow(x,2) + -3 * x + 7 + 4 * np.random.rand(1,150)
xx = open("data_150.txt","w")
xx.write("150")
for i in range(150):
xx.write("\n")
xx.write(str(round(x[0][i],5)))
xx.write("\t")
xx.write(str(round(y[0][i],5)))
xx.close()
import pandas as pd
xxx = open("3-dgree-regression_results.txt", "r")
data = xxx.read()
data = data.split("\n")
n_data = []
for i in data:
n_data.append(i.split("\t"))
data = n_data[:-1]
X = []
Y = []
for i in data:
X.append(float(i[0]))
Y.append(float(i[1]))
plt.title("data for regression")
plt.plot(X,Y,"r-",label="recovered")
plt.plot(x,y, "b.")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()
Raw
polynomialRegiression.cpp
// https://himbopsa.tistory.com/13
#include <iostream>
#include <fstream>
#include <Eigen/Dense>
#include <cmath>
#include <Eigen/LU>
#include <chrono>
using namespace Eigen;
using namespace std;
int main(int, char**) {
float N; // the number of dates
ifstream file("/Users/newpolaris/Projects/temp/temp/data_150.txt");
if (!file) {
cout << "file path error";
return -1;
}
file >> N;
float *data_x = new float[N];
float *data_y = new float[N];
for (int i = 0; i < N; i++) {
file >> data_x[i] >> data_y[i];
cout << data_x[i];
cout << " " << data_y[i] << endl;
}
float Nx = 0, Nx2 = 0, Nx3 = 0, Nx4 = 0,
Nx5 = 0, Nx6 = 0, y = 0, yx = 0,
yx2 = 0, yx3 = 0;
Matrix4d X, XI;
Vector4d Y, O, O2, O3, O4, O5;
for (int i = 0; i < N; i++) {
Nx += pow(data_x[i], 1.f);
Nx2 += pow(data_x[i], 2.f);
Nx3 += pow(data_x[i], 3.f);
Nx4 += pow(data_x[i], 4.f);
Nx5 += pow(data_x[i], 5.f);
Nx6 += pow(data_x[i], 6.f);
y += data_y[i];
yx += data_y[i] * pow(data_x[i], 1.f);
yx2 += data_y[i] * pow(data_x[i], 2.f);
yx3 += data_y[i] * pow(data_x[i], 3.f);
}
X(0, 0) = Nx3; X(0, 1) = Nx2; X(0, 2) = Nx; X(0, 3) = N;
X(1, 0) = Nx4; X(1, 1) = Nx3; X(1, 2) = Nx2; X(1, 3) = Nx;
X(2, 0) = Nx5; X(2, 1) = Nx4; X(2, 2) = Nx3; X(2, 3) = Nx2;
X(3, 0) = Nx6; X(3, 1) = Nx5; X(3, 2) = Nx4; X(3, 3) = Nx3;
Y[0] = y; Y[1] = yx; Y[2] = yx2; Y[3] = yx3;
// cout << X << endl;
// cout << Y << endl;
MatrixXd Xd(X);
chrono::steady_clock clock;
auto t0 = clock.now();
// https://eigen.tuxfamily.org/dox-devel/group__LeastSquares.html
O = X.inverse() * Y;
auto t1 = clock.now();
O2 = (X.transpose() * X).ldlt().solve(X.transpose() * Y);
auto t2 = clock.now();
O3 = X.colPivHouseholderQr().solve(Y);
auto t3 = clock.now();
O4 = Xd.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y);
auto t4 = clock.now();
cout << O[0] << "x3" << " + " << O[1] << "x2"
<< O[2] << "x1" << " + " << O[3] << endl;
cout << O2[0] << "x3" << " + " << O2[1] << "x2"
<< O2[2] << "x1" << " + " << O2[3] << endl;
cout << O3[0] << "x3" << " + " << O3[1] << "x2"
<< O3[2] << "x1" << " + " << O3[3] << endl;
cout << O4[0] << "x3" << " + " << O4[1] << "x2"
<< O4[2] << "x1" << " + " << O4[3] << endl;
chrono::duration<float, std::milli> t[4];
t[0] = t1 - t0;
t[1] = t2 - t1;
t[2] = t3 - t2;
t[3] = t4 - t3;
file.close();
for (int i = 0; i < 4; i++)
cout << "t[" << i << "] : " << t[i].count() << "\n";
ofstream result("3-dgree-regression_results.txt");
cout << "x : ";
const float step = 0.0033;
for (int i = 0; i < N; i++) {
const float x = round(step * i * 10000) / 1000;
const float x2 = x*x, x3 = x*x*x;
cout << x << " ";
result << x << "\t" << O[0] * x3 + O[1] * x2 + O[2] * x + O[3];
result << endl;
}
cout << endl;
result.close();
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment