In [1]:
from numpy import *

In [3]:
n = 3
x = random.rand(n)
y = 3*x - 5 + 0.2*random.randn(n)
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x,y,'bo');

In [7]:
def E1(alpha,beta):
global x,y
e = alpha*x + beta - y  # array of pointwise errors
return abs(e).sum()

def E2(alpha,beta):
global x,y
e = alpha*x + beta - y  # array of pointwise errors
return (e**2).sum()

# try a bunch of points (alpha,beta) in some box that we think contains the minimizer ...
p = 500000
alphas =  2 + 2*random.rand(p)  # in interval (2,4)
betas  = -6 + 2*random.rand(p)  # in interval (-6,-4)
# try each (alpha,beta) in turn and keep the best

Emin = inf  # best (least) we've seen so far
for alpha,beta in zip(alphas,betas):
E = E1(alpha,beta)
if E < Emin:
Emin = E
alpha1,beta1 = alpha,beta  # best choice seen so far
print(alpha1,beta1)
plt.plot(x,y,'bo')
xx = array([0,1])
yy = alpha1*xx + beta1
plt.plot(xx,yy,'r')
plt.savefig('temp.svg')
import os
os.system('firefox temp.svg')

3.16480410932 -5.10688867997

Out[7]:
0

In [9]:
7+linspace(0,10,11)

Out[9]:
array([  7.,   8.,   9.,  10.,  11.,  12.,  13.,  14.,  15.,  16.,  17.])
In [ ]:
Examples

In [10]:
a = ones((3,4))
b = ones((3,4))
a+b

Out[10]:
array([[ 2.,  2.,  2.,  2.],
[ 2.,  2.,  2.,  2.],
[ 2.,  2.,  2.,  2.]])
In [12]:
a = random.rand(3,4)
b = ones((3,1))
a+b

Out[12]:
array([[ 1.63473869,  1.32671351,  1.40863303,  1.45192962],
[ 1.88646835,  1.09996372,  1.83734529,  1.11091262],
[ 1.53776573,  1.90239007,  1.71366127,  1.74806915]])
In [13]:
a = random.rand(3,4,2)
b = ones((3,1,2))
a+b

Out[13]:
array([[[ 1.73647591,  1.45934946],
[ 1.63406906,  1.33354595],
[ 1.62710045,  1.25239646],
[ 1.43293443,  1.95587124]],

[[ 1.00304175,  1.30599109],
[ 1.59972272,  1.41963777],
[ 1.68986117,  1.4773109 ],
[ 1.48119428,  1.85920368]],

[[ 1.14412625,  1.95004026],
[ 1.91910986,  1.29182713],
[ 1.31333243,  1.61809557],
[ 1.90693232,  1.70454338]]])
In [19]:
alpha = array([1,3,5])
x = array([1,10,100,1000])
x2 = x.reshape( (len(x),1) )
print(alpha)
print(x)
print(x2)
foo = (alpha*x2).T
print(foo)

[1 3 5]
[   1   10  100 1000]
[[   1]
[  10]
[ 100]
[1000]]
[[   1   10  100 1000]
[   3   30  300 3000]
[   5   50  500 5000]]

In [17]:
foo.sum()

Out[17]:
9999
In [20]:
foo.sum(axis=1)  # sum over 2nd index, i.e. across columns

Out[20]:
array([1111, 3333, 5555])
In [50]:
random.seed(6436)
n = 3
x = random.rand(n)
y = 3*x - 5 + 0.1*random.rand(n)

p = 50000
alphas =  2 + 2*random.rand(p)  # in interval (2,4)
betas  = -6 + 2*random.rand(p)  # in interval (-6,-4)
# try each (alpha,beta) in turn and keep the best

alphas2 = alphas.reshape((p,1))
ax = alphas2*x
axT = ax.T
axTb = axT + betas
#print(axTb, axTb.shape)
axTbT = axTb.T
e = axTbT - y
E = abs(e).sum(axis=1)
#print(E)
i = E.argmin()
#print(i)
alpha1 = alphas[i]
beta1 = betas[i]
plt.plot(x,y,'bo')
xx = array([0,1])
yy = alpha1*xx + beta1
plt.plot(xx,yy,'r')
plt.savefig('temp.svg')
import os
os.system('firefox temp.svg')

Out[50]:
0
In [51]:
plt.scatter(alphas,betas,c=E,edgecolor='none')
plt.colorbar()
plt.xlabel('alpha')
plt.ylabel('beta');