  We extend to large scales a method proposed in previous work that
reconstructs non-parametrically the primordial power spectrum from cosmic
microwave background data at high resolution. The improvement is necessary to
account for the non-gaussianity of the Wilkinson Microwave Anisotropy Probe
(WMAP) likelihood due primarily to cosmic variance. We assume the concordance
LambdaCDM cosmology, utilise a smoothing prior and perform Monte Carlo
simulations around an initial power spectrum that is scale-free and with
spectral index n_s=0.97, very close to the concordance spectrum. The horizon
scale for the model we are considering corresponds to the wavenumber
k_h=4.52X10^{-4} Mpc^{-1}. We find some evidence for the presence of features
and we quantify the probabilities of exceeding the observed deviations in WMAP
data with respect to the fiducial models. We detect the following marginal
departures from a scale-free (spectral index n_s=0.97) initial spectrum: a
cut-off at 0.0001<k<0.001 Mpc^{-1} at 79.5% (92%), a dip at 0.001<k<0.003
Mpc^{-1} at 87.2% (98%) and a bump at 0.003<k<0.004 Mpc^{-1} at 90.3% (55.5%)
confidence level.
  These frequentist confidence levels are calculated by integrating over the
distribution of the Monte Carlo reconstructions built around the fiducial
models. The frequentist analysis finds the low k cutoff of the estimated power
spectrum to be about 2.5 sigma away from the n_s=0.97 model, while in the
Bayesian analysis the model is about 1.5 sigma away from the estimated
spectrum. (The sigma's are different for the two different methods.)
