Given Position Measurements, How to Estimate Velocity and Acceleration
One approach would be to cast the problem as least-squares smoothing. The idea is to locally fit a polynomial with a moving window, then evaluate the derivative of the polynomial. This answer about Savitzky-Golay filtering has some theoretical background on how it works for nonuniform sampling. 

In this case, code is probably more illuminating as to the benefits/limitations of the technique. The following numpy script will calculate the velocity and acceleration of a given position signal based on two parameters: 1) the size of the smoothing window, and 2) the order of the local polynomial approximation.

# Example Usage:
# python sg.py position.dat 7 2

import math
import sys

import numpy as np
import numpy.linalg
import pylab as py

def sg_filter(x, m, k=0):
    """
    x = Vector of sample times
    m = Order of the smoothing polynomial
    k = Which derivative
    """
    mid = len(x) / 2        
    a = x - x[mid]
    expa = lambda x: map(lambda i: i**x, a)    
    A = np.r_[map(expa, range(0,m+1))].transpose()
    Ai = np.linalg.pinv(A)

    return Ai[k]

def smooth(x, y, size=5, order=2, deriv=0):

    if deriv &gt; order:
        raise Exception, "deriv must be &lt;= order"

    n = len(x)
    m = size

    result = np.zeros(n)

    for i in xrange(m, n-m):
        start, end = i - m, i + m + 1
        f = sg_filter(x[start:end], order, deriv)
        result[i] = np.dot(f, y[start:end])

    if deriv &gt; 1:
        result *= math.factorial(deriv)

    return result

def plot(t, plots):
    n = len(plots)

    for i in range(0,n):
        label, data = plots[i]

        plt = py.subplot(n, 1, i+1)
        plt.tick_params(labelsize=8)
        py.grid()
        py.xlim([t[0], t[-1]])
        py.ylabel(label)

        py.plot(t, data, 'k-')

    py.xlabel("Time")

def create_figure(size, order):
    fig = py.figure(figsize=(8,6))
    nth = 'th'
    if order &lt; 4:
        nth = ['st','nd','rd','th'][order-1]

    title = "%s point smoothing" % size
    title += ", %d%s degree polynomial" % (order, nth)

    fig.text(.5, .92, title,
             horizontalalignment='center')

def load(name):
    f = open(name)    
    dat = [map(float, x.split(' ')) for x in f]
    f.close()

    xs = [x[0] for x in dat]
    ys = [x[1] for x in dat]

    return np.array(xs), np.array(ys)

def plot_results(data, size, order):
    t, pos = load(data)
    params = (t, pos, size, order)

    plots = [
        ["Position",     pos],
        ["Velocity",     smooth(*params, deriv=1)],
        ["Acceleration", smooth(*params, deriv=2)]
    ]

    create_figure(size, order)
    plot(t, plots)

if __name__ == '__main__':
    data = sys.argv[1]
    size = int(sys.argv[2])
    order = int(sys.argv[3])

    plot_results(data, size, order)
    py.show()


Here are some example plots (using the data you provided) for various parameters.







Note how the piecewise constant nature of the acceleration becomes less obvious as the window size increases, but can be recovered to some extent by using higher-order polynomials. Of course, other options involve applying a first derivative filter twice (possibly of different orders). Another thing that should be obvious is how this type of Savitzky-Golay filtering, since it uses the midpoint of the window, truncates the ends of the smoothed data more and more as the window size increases. There are various ways to solve that problem, but one of the better ones is described in the following paper:


  P.A. Gorry, General least-squares smoothing and differentiation by the
  convolution (Savitzky–Golay) method, Anal. Chem. 62 (1990) 570–573. (google)


Another paper by the same author describes a more efficient way to smooth nonuniform data than the straightforward method in the example code:


  P.A. Gorry, General least-squares smoothing and differentiation of
  nonuniformly spaced data by the convolution method, Anal. Chem. 63
  (1991) 534–536. (google)


Finally, one more paper worth reading in this area is by Persson and Strang:


  P.O. Persson, G. Strang, Smoothing by Savitzky–Golay and Legendre
  Filters, Comm. Comp. Finance 13 (2003) 301–316. (pdf link)


It contains much more background theory, and concentrates on error analysis for choosing a window size.