We study the sample variance of the matter power spectrum for the standard Lambda Cold Dark Matter universe. We use a total of 5000 cosmological N-body simulations to study in detail the distribution of best-fit cosmological parameters and the baryon acoustic peak positions. The obtained distribution is compared with the results from the Fisher matrix analysis with and without including non-Gaussian errors. For the Fisher matrix analysis, we compute the derivatives of the matter power spectrum with respect to cosmological parameters using directly full nonlinear simulations. We show that the non-Gaussian errors increase the unmarginalized errors by up to a factor 5 for k_{max}=0.4h/Mpc if there is only one free parameter provided other parameters are well determined by external information. On the other hand, for multi-parameter fitting, the impact of the non-Gaussian errors is significantly mitigated due to severe parameter degeneracies in the power spectrum. The distribution of the acoustic peak positions is well described by a Gaussian distribution, with its width being consistent with the statistical interval predicted from the Fisher matrix. We also examine systematic bias in the best-fit parameter due to the non-Gaussian errors. The bias is found to be smaller than the 1 sigma statistical error for both the cosmological parameters and the acoustic scale positions. ; Comment: 12 pages, 10 figures, accepted for publication in ApJ, minor changes