In [21]:

    

import numpy as np
from numpy import ones, hstack, matrix, divide, multiply
from numpy.random import random as nprandom

M = hstack((ones((100, 1)), (100 * nprandom((100, 4)))))
X = M[:, 1:2]
Y = M[:, 2:3]
Z = M[:, 3:4]

# F = 3x²/y - 17xy + 8/z + 0w - 13

X2_Y = divide(multiply(X, X), Y)
XY = multiply(X, Y)
_Z = 1 / Z

M2 = hstack((M, X2_Y, XY, _Z))
C = np.matrix('-13; 0; 0; 0; 0; 3; -17; 8')

F = M2 * C  # correct values of F

TETA = np.linalg.inv(M2.transpose().dot(M2)).dot(M2.transpose()).dot(F)
# TETA = inverse(M2' * M2) * M2' * F; # normal equation

print TETA - C
print 'Holy chupacabra!'


    

    




[[  5.15319698e-09]
 [ -8.12576673e-11]
 [ -1.90993887e-11]
 [  7.95807864e-13]
 [  4.86188867e-12]
 [ -1.33226763e-14]
 [  2.62900812e-13]
 [  5.20230969e-10]]
Holy chupacabra!