gistfile1.txt

eta = np.array([1, 0.01]);
theta0 = 0.5;
N = 30;
import numpy as np
xl = linspace(-1,1,N)
xs = np.zeros((2, N*N));

for x1 in xrange(0,N):
    for x2 in xrange(0, N):
        j = x2*N + x1;
        xs[:, j] = np.array([x1, x2]);

def ard_kernel(x1, x2):
    sum = eta[0] * (x1[0] - x2[0])**2 + eta[1] * (x1[1] - x1[1])**2
    retval = theta0 * exp(-1/2.0*sum)
    return retval

from scipy import linalg as slinalg
def big_rand(mu, cov):
    N = len(mu);
    L = slinalg.cholesky(cov, lower=True)
    z = randn(N,1)
    R = mu + L * z
    return R

def outer_f( x, y, func):
     xl = len(x[0])
     yl = len(y[0])
     q = np.zeros((xl,yl))
     for xi in xrange(xl):
        for yi in xrange(yl):
            q[xi, yi] = func( x[:,xi], y[:, yi] )
     return q

cov = outer_f(xs, xs, ard_kernel)

if bm < 0:
    beta = -bm
else:
    beta = 1.0e-14
cov += np.eye(N*N) * beta;

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = Axes3D(fig)
R = big_rand(zeros((N*N, 1)), cov);
xx, yy = np.meshgrid(xl, xl);
# z = np.sin(xx) + np.cos(yy)
ax.plot_wireframe(xx, yy, z)
plt.show()