Skip to content

Instantly share code, notes, and snippets.

@matsuken92
Last active August 29, 2015 14:15
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save matsuken92/1fbd89326c26a15e0442 to your computer and use it in GitHub Desktop.
Save matsuken92/1fbd89326c26a15e0442 to your computer and use it in GitHub Desktop.
pythonを使った回帰分析の概念の解説 その2 ref: http://qiita.com/kenmatsu4/items/1e37da1d55292035d985
S(\alpha, \beta) =
\left( \sum_i^n x_i^2 \right) \alpha^2 + n\beta^2
+ 2 \left( \sum_i^n x_i \right)\alpha \beta
- 2 \left( \sum_i^n x_i y_i \right)\alpha
- 2 \left( \sum_i^n y_i \right)\beta
+ \sum_i^n y_i^2
S(\beta) = n\beta^2
+ 2 \left( \sum_i^n (x_i\alpha - y_i) \right) \beta
+ \alpha^2\sum_i^n x_i^2 - 2\alpha \sum_i^n x_iy_i + \sum_i^n y_i^2
\frac{\partial S}{\partial \alpha} = 0,
\frac{\partial S}{\partial \beta } = 0,
\frac{\partial S}{\partial \alpha} = 2\left(\sum_i^n x_i^2 \right) \alpha + 2\left( \sum_i^n x_i \right) \beta - 2\sum_i^n x_i y_i = 0
\frac{\partial S}{\partial \beta} = 2n\beta + 2\left( \sum_i^n x_i\right) \alpha - 2\sum_i^ny_i = 0
\frac{1}{2}\frac{\partial S}{\partial \alpha} = 34288 \alpha + 1240 \beta - 18884 = 0
\frac{1}{2}\frac{\partial S}{\partial \beta} = 50 \beta + 1240 \alpha - 655 = 0
from sympy import *
a, b = symbols('a b')
init_printing()
# 34288α + 1240β − 18884 = 0
# 50β + 1240α − 655 = 0
solve([34288 * a + 1240 * b - 18884, 50* b + 1240 * a - 655], [a, b])
a --> 0.746606334842
b --> -5.41583710407
\alpha = 0.746606334842
\beta = -5.41583710407
import numpy as np
import matplotlib.pyplot as plt
data= np.loadtxt('cars.csv',delimiter=',',skiprows=1)
data[:,1] = map(lambda x: x * 1.61, data[:,1]) # mph から km/h に変換
data[:,2] = map(lambda y: y * 0.3048, data[:,2]) # ft から m に変換
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
ax.set_xlim(0,50)
ax.set_title("Stopping Distances of Cars with estimated regression line")
ax.set_xlabel("speed(km/h)")
ax.set_ylabel("distance(m)")
plt.scatter(data[:,1],data[:,2])
x = np.linspace(0,50,50)
y = 0.746606334842 * x -5.41583710407
plt.plot(x,y)
\left( \sum_i^n x_i^2 \right) =34288
\left( \sum_i^n y_i^2 \right)=11604
\left( \sum_i^n x_iy_i \right)=18884
\left( \sum_i^n x_i \right)=1240
\left( \sum_i^n y_i \right)=655
S(\alpha, \beta) =
34288 \alpha^2 + 50\beta^2
+ 2480\alpha \beta
- 37768\alpha
- 1310 \beta
+ 11604
from mpl_toolkits.mplot3d.axes3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
# set field
X = np.linspace(0.2, 1.3, 100)
Y = np.linspace(-25, 15, 100)
# set data
#sum(x**2)
sum_x_2 = 34288.2988
#sum(y**2)
sum_y_2 = 11603.8684051
#dot(x,y)
sum_xy = 18884.194896
#sum(x)
sum_x = 1239.7
#sum(y)
sum_y = 655.0152
X, Y = np.meshgrid(X, Y)
#S(α,β)=34288α^2 + 50β^2 + 2480αβ − 37768α − 1310β + 11604
S = (sum_x_2 * (X**2)) + (50 * (Y**2)) + (2 * sum_x * X * Y) + (-2 * sum_xy * X) + (-2 * sum_y * Y) + sum_y_2
# prepare plot
fig = plt.figure(figsize=(18,6))
ax = fig.add_subplot(121, projection='3d', azim=60)
ax.set_xlabel("alpha")
ax.set_ylabel("beta")
ax.set_zlabel("S")
# draw 3D graph
surf = ax.plot_surface(X, Y, S, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
# draw contour
ax = fig.add_subplot(122)
plt.contour(X,Y,S,50)
ax.set_xlabel("alpha")
ax.set_ylabel("beta")
plt.show()
S(\alpha) = \left( \sum_i^n x_i^2 \right) \alpha^2
+ 2\left( \sum_i^n (x_i\beta - x_i y_i ) \right) \alpha
+ n\beta^2 - 2\beta\sum_i^n y_i + \sum_i^n y_i^2
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment