Why can't I reproduce Hinton's results ?

I am trying to reproduce Hinton's results in his paper Improving neural networks by preventing co-adaptation of feature detectors. With a two hidden layers architecture, with 50% dropout in the hidden layers and 800 neurons in each, he achieve a 1.3% error rate on the MNIST database.

I tried to do the same using Pylearn2. I monitor using cross-validation, and stop when the score haven't been increased for the past 100 iterations. It run for 683 iteration. The final error rate on the test database is 2.02%.

I just found an old post of Ian Goodfellow about it, but there is no indication that he succeded.

Here is my code :

from __future__ import division
import os

import numpy as np

from pylearn2.train import Train
from pylearn2.datasets.mnist import MNIST
from pylearn2.models import softmax_regression, mlp
from pylearn2.training_algorithms import bgd, sgd
from pylearn2.termination_criteria import MonitorBased
from pylearn2.train_extensions import best_params
from pylearn2.utils import serial
from pylearn2.costs.mlp.dropout import Dropout
from theano import function
from theano import tensor as T

h0 = mlp.Softplus(layer_name='h0', dim=800, sparse_init=40)
h1 = mlp.Softplus(layer_name='h1', dim=800, sparse_init=40)
#h2 = mlp.Softplus(layer_name='h2', dim=50, sparse_init=15)
ylayer = mlp.Softmax(layer_name='y', n_classes=10, irange=0)
layers = [h0, h1, ylayer]

model = mlp.MLP(layers, nvis=784)
train = MNIST('train', one_hot=1, start=0, stop=50000)
valid = MNIST('train', one_hot=1, start=50000, stop=60000)
test = MNIST('test', one_hot=1, start=0, stop=10000)

monitoring = dict(valid=valid)
termination = MonitorBased(channel_name="valid_y_misclass", N=100)
extensions = [best_params.MonitorBasedSaveBest(channel_name="valid_y_misclass",
save_path="train_best.pkl")]

algorithm = sgd.SGD(0.1, batch_size=100, cost=Dropout(),
                    monitoring_dataset = monitoring, termination_criterion = termination)

save_path = "train_best.pkl"
if os.path.exists(save_path):
    model = serial.load(save_path)
else:
    print 'Running training'
    train_job = Train(train, model, algorithm, extensions=extensions, save_path="train.pkl", save_freq=1)
    train_job.main_loop()

X = model.get_input_space().make_batch_theano()
Y = model.fprop(X)

y = T.argmax(Y, axis=1)
f = function([X], y)
yhat = f(test.X)

y = np.where(test.get_targets())[1]
print 'accuracy', (y==yhat).sum() / y.size