You are an expert at summarizing long articles. Proceed to summarize the following text: 
millisecond pulsars ( msps ) have for some time been known to exhibit exceptional rotational stability , with decade long observations providing timing measurements with accuracies similar to atomic clocks ( e.g. ) . such stability lends itself well to the pursuit of a wide range of scientific goals , e.g. observations of the pulsar psr b1913 + 16 showed a loss of energy at a rate consistent with that predicted for gravitational waves @xcite , whilst the double pulsar system psr j0737 - 3039a / b has provided precise measurements of several ` post keplerian ' parameters allowing for additional stringent tests of general relativity @xcite .    for a detailed review of pulsar timing refer to e.g. @xcite . in brief , the arrival times of pulses ( toas ) for a particular pulsar will be recorded by an observatory in a series of discrete observations over a period of time . these arrival times must all be transformed into a common frame of reference , the solar system barycenter , in order to correct for the motion of the earth .    a model for the pulsar can then be fitted to the toas ; this characterises the properties of the pulsar s orbital motion , as well as its timing properties such as its orbital frequency and spin down . this is most commonly carried out using the tempo2 pulsar - timing packages @xcite , or more recently , the bayesian pulsar timing package temponest @xcite .    when performing this fitting process , both tempo2 and temponest assume purely gaussian statistics in the properties of the uncorrelated noise . in realistic datasets , however , this assumption is not necessarily correct . if the underlying probability density function ( pdf ) for the noise is not gaussian , for example , if there is an excess of outliers relative to a gaussian distribution , modifiers to the toa error bars that scale their size are used to find the best approximation to a gaussian distribution . this can be performed using a single modifier for a given receiving system determined across an entire dataset , or as in the ` fixdata ' plugin for tempo2 @xcite , where the modifier is determined separately for a series of short time lags . while the latter of these two approaches can better account for a non - gaussian distribution in the uncorrelated noise , it does so at the expense of a potentially large number of additional free parameters , and ultimately does not address the core issue , that the underlying distribution is not gaussian . both approaches then have the direct consequence of decreasing the precision with which one can estimate the timing parameters , and any other signals of interest , such as intrinsic red spin noise due to rotational irregularities in the neutron star @xcite or correlated noise due to a stochastic gravitational wave background ( gwb ) generated by , for example , coalescing black holes ( e.g. @xcite ) . indeed , currently all published limits on the signals induced by a gwb have been obtained under the assumption that the statistics of the toa errors are gaussian ( see e.g. @xcite ) .    in this paper we introduce a method of performing a robust bayesian analysis of non - gaussianity present in pulsar timing data , simultaneously with the pulsar timing model , and additional stochastic parameters such as those describing the red noise , and dispersion measure variations . the parameters used to define the presence of non - gaussianity are zero for gaussian processes , giving a simple method of defining the strength of non - gaussian behaviour . in section [ section : bayes ] we will describe the basic principles of our bayesian approach to data analysis , giving a brief overview of how it may be used to perform model selection , and introduce multinest . in sections [ section : nongausslike ] and [ section : toy ] we introduce the non - gaussian likelihood we will use in our pulsar timing analysis , and apply it to a simple toy problem . in section [ section : pulsarnongaussian ] we then extend this likelihood to the subject of pulsar timing , and apply it to both simulated and real data in sections [ section : pulsarsims ] and [ section : realdata ] respectively , before finally offering some concluding remarks in section [ section : conclusion ] . this research is the result of the common effort to directly detect gravitational waves using pulsar timing , known as the european pulsar timing array ( epta ) @xcite . given a set of data @xmath1 , bayesian inference provides a consistent approach to the estimation of a set of parameters @xmath2 in a model or hypothesis @xmath3 . in particular , bayes theorem states that :    @xmath4    where @xmath5 is the posterior probability distribution of the parameters , @xmath6 is the likelihood , @xmath7 is the prior probability distribution , and @xmath8 is the bayesian evidence . since the evidence is independent of the parameters @xmath2 it is typically ignored when one is only interested in performing parameter estimation . in this case inferences are obtained by taking samples from the ( unnormalised ) posterior using , for example , standard markov chain monte carlo ( mcmc ) sampling methods .    for model selection , however , the evidence is key , and is defined simply as the factor required to normalise the posterior over @xmath2 :    @xmath9    where @xmath10 is the dimensionality of the parameter space .    as the evidence is just the average of the likelihood over the prior , it will be larger for a simpler model with a compact parameter space if more of that parameter space is likely . more complex models where large areas of parameter space have low likelihood values will have a smaller evidence even if the likelihood function is very highly peaked , unless they are significantly better at explaining the data . thus , the evidence automatically implements occam s razor . the question of model selection between two models @xmath11 and @xmath12 can be answered via the model selection ratio @xmath13 , commonly referred to as the ` bayes factor ' :    @xmath14    where @xmath15 is the a priori probability ratio for the two models , which in this work we will set to unity but occasionally requires further consideration . the bayes factor then allows us to obtain the probability of one model compared the other simply as :    @xmath16    in practice when performing bayesian analysis we do not work with the likelihood , but the log likelihood . in this case the quantity of interest is the log bayes factor , which is simply the difference in the log evidence for the two models . for example , a difference in the log evidence of 3 for two competing models gives a bayes factor of @xmath17 , which in turn gives a probability of @xmath18 . we use the difference in the log evidence in sections [ section : pulsarsims ] and [ section : realdata ] to perform model selection between our gaussian and non - gaussian models . while many techniques exist for calculating the evidence , such as thermodynamic integration @xcite , it remains a challenging task both numerically and computationally , with evidence evaluation at least an order - of - magnitude more costly than parameter estimation . nested sampling @xcite is an approach designed to make the calculation of the evidence more efficient , and also produces posterior inferences as a by - product . the multinest algorithm ( @xcite ) builds upon this nested sampling framework , and provides an efficient means of sampling from posteriors that may contain multiple modes and/or large ( curving ) degeneracies , and also calculates the evidence . since its release multinest has been used successfully in a wide range of astrophysical problems , including inferring the properties of a potential stochastic gravitational wave background in pulsar timing array data @xcite , and is also used in the bayesian pulsar timing package temponest . this technique has greatly reduced the computational cost of bayesian parameter estimation and model selection , and is employed in this paper . in this section we will outline the method adopted for including non - gaussian behaviour in our analysis . we use the approach developed in @xcite , which is based on the energy eigenmode wavefunctions of a simple harmonic oscillator . we will describe this in brief below in order to aid future discussion . we begin by considering our data , the vector @xmath19 of length @xmath20 , as the sum of some signal @xmath21 and noise @xmath22 such that : @xmath23 we can then construct the likelihood that the residuals after subtracting our model signal from the data follows an uncorrelated gaussian distribution of width @xmath24 as :    @xmath25,\ ] ]    with @xmath26 the diagonal noise covariance matrix for the residuals , such that @xmath27 , and @xmath28 the determinant of @xmath26 . we now extend this to the general case in order to allow for non - gaussian distributions by modelling our pdf as the sum of a set of gaussians , modified by hermite polynomials @xmath29 ( see e.g. @xcite for previous uses of hermite polynomials in describing departures from gaussianity ) , defined as : @xmath30    therefore , for a general random variable @xmath31 the pdf for fluctuations in @xmath31 can be written : @xmath32\left|\sum_{n=0}^{\infty}\alpha_nc_nh_n\left(\frac{x}{\sqrt{2}\sigma}\right)\right|^2\ ] ] with @xmath33 free parameters that describe the relative contributions of each term to the sum , and @xmath34 is a normalization factor . equation [ eq : prx ] forms a complete set of pdfs , normalised such that :    @xmath35c_nh_n\left(\frac{x}{\sqrt{2}\sigma}\right)c_mh_m\left(\frac{x}{\sqrt{2}\sigma}\right ) = \delta_{mn},\ ] ]    with @xmath36 the kronecker delta , where the ground state , @xmath11 , reproduces a standard gaussian pdf , and any non - gaussianity in the distribution of @xmath31 will be reflected in non - zero values for the coefficients @xmath33 associated with higher order states . the only constraint we must place on the values of the amplitudes @xmath37 is : @xmath38 with @xmath39 the maximum number of coefficients to be included in the model for the pdf . this is performed most simply by setting : @xmath40    we can therefore rewrite eq . [ eq : gausslike ] in this more general form as : @xmath41\nonumber \\ & \times&\prod_{i=1}^{n_d}\left|\sum_{n=0}^{n_{\mathrm{max}}}\alpha_nc_nh_n\left(\frac{d_i - s_i}{\sqrt{2}\sigma}\right)\right|^2.\end{aligned}\ ] ] the advantage of this method is that one may use a finite set of non - zero @xmath33 to model the non - gaussianity , without mathematical inconsistency . any truncation of the series still yields a proper distribution , in contrast to the more commonly used edgeworth expansion ( e.g. @xcite ) . before applying the formalism described in section [ section : nongausslike ] to the practice of pulsar timing , we first demonstrate its use in a toy problem . here our data vector @xmath19 contains 10000 points drawn from a non - gaussian distribution obtained using eq . [ eq : prx ] , with parameters listed in table [ table : toyparams ] . .parameters used to generate non - gaussian noise in a simple toy problem . [ cols="^,^,^ " , ]     [ table : logevidence ]    table [ table : logevidence ] lists the log evidence values for different sets of non - gaussian coefficients , normalised such that the log evidence for no additional coefficients ( i.e. assuming gaussian statistics ) is 0 . we see that there is a significant increase in the log evidence ( @xmath42 39 ) when including even just two coefficients , indicating definitive support for their inclusion in the model . as the number increases the rise in evidence increases , reaching a maximum with 4 included coefficients . given the timing model , red noise and dispersion measure variation solutions that were subtracted from the data were obtained from a gaussian analysis , we will however still include coefficients up to and including @xmath43 in the full analysis . given the large dimensionality of the problem this analysis can not be carried out using multinest . as such we make use of the guided hamiltonian sampler used previously in pulsar timing analysis in @xcite . this sampler makes use of both gradient information in the likelihood , and also the hessian in order to efficiently sample from large parameter spaces . table [ table:0437 ] lists the timing model parameter estimates and their nominal standard deviations for both the gaussian and non - gaussian analysis . in all cases we find the parameter estimates and their uncertainties to be consistent between both methods . in fig . [ figure:0437noise ] ( top ) we show the one and two - dimensional marginalised posterior distributions for the red noise and dispersion measure variation power law amplitudes and spectral indicies for the non - gaussian ( left ) and gaussian ( right ) analysis . both are also extremely consistent with one another , however when overlaying the two sets of 1-dimensional posterior distributions for each of the 4 parameters separately ( bottom 4 panels ) some differences become apparent between the non - gaussian ( blue dashed lines ) and gaussian ( red solid lines ) analysis . in particular the dispersion measure variation power law parameter estimates show a slight shift towards higher amplitudes and shallower spectral indices in the non - gaussian case .    despite these similarities in the timing and stochastic parameter estimates between the gaussian and non - gaussian analysis , fig . [ figure:0437pdf ] indicates a definitive detection of non - gaussianity in the dataset , in agreement with the difference in the log evidence for the noise only analysis . in the top plot we show the one and two - dimensional marginalised posterior distributions for the 5 non - gaussian coefficients fit in the analysis of j0437@xmath04715 . vertical lines are included at 0 where visible in the plots , however , except for @xmath44 all the coefficients are inconsistent with this value . in the bottom plot we then show the set of equally weighted pdfs obtained from the non - gaussian analysis ( black lines ) setting @xmath45 . in addition we over plot the mean of the distribution ( red line ) and a unit gaussian ( blue line ) all of which have been normalised to have a sum of 1 . the difference between the gaussian and non - gaussian pdfs is clear , with a larger probability for both small ( @xmath46 ) and larger ( @xmath47 ) deviations than given by the gaussian pdf . that such a significant detection of non - gaussianity does not lead to larger changes in the parameter estimates can potentially be attributed to a frequency dependence on the significance of the @xmath48 parameters . in @xcite the 10 cm j0437@xmath04715 data was found to be describable through gaussian statistics alone . this would suggest that the non - gaussianity we detect exists primarily at low frequencies . in figure [ figure:0437normres ] we show the normalised residuals from fig . [ figure:0437noise ] separated into its 10 cm , 20 cm and 50 cm components , along with histograms for each wavelength . here the increase in non - gaussian behaviour can clearly be seen as the wavelength increases . given the lowest frequencies have the greatest degree of non - gaussianity it is less surprising that there is little impact on the timing or red spin noise parameters , as the low frequency data contributes the least to these parts of the model . the low frequencies do , however , contribute greatly to the constraints on dispersion measure variations , and it is here we see the greatest difference between the gaussian and non - gaussian models . lcc + pulsar name & j0437@xmath04715 + mjd range & 50191.055619.2 + data span ( yr ) & 14.86 + number of toas & 5052 +   + model parameter & non  gaussian & gaussian + right ascension , @xmath48 ( rad ) & 1.20979650940(10 ) & 1.20979650943(11 ) + declination , @xmath49 ( rad ) & @xmath00.82471224153(8 ) & @xmath00.82471224154(8 ) + pulse frequency , @xmath50 ( s@xmath51 ) & 173.6879458121850(3 ) & 173.6879458121849(4 ) + first derivative of pulse frequency , @xmath52 ( s@xmath53 ) & @xmath01.728365(4)@xmath54 & @xmath01.728365(4)@xmath54 + dispersion measure , dm ( @xmath55pc ) & 2.64462(11 ) & 2.64461(11 ) + first derivative of dispersion measure , @xmath56 ( @xmath55pcyr@xmath51 ) & @xmath06(6)@xmath57 & @xmath07(7)@xmath57 + dm2 ( @xmath55 pc yr@xmath53 ) & @xmath01(2)@xmath58 & @xmath01(2)@xmath58 + proper motion in right ascension , @xmath59 ( masyr@xmath51 ) & 121.439(3 ) & 121.441(3 ) + proper motion in declination , @xmath60 ( masyr@xmath51 ) & @xmath071.474(3 ) & @xmath071.474(3 ) + parallax , @xmath61 ( mas ) & 6.4(2 ) & 6.3(2 ) + orbital period , @xmath62 ( d ) & 5.7410462(3 ) & 5.7410461(3 ) + epoch of periastron , @xmath63 ( mjd ) & 54530.1722(3 ) & 54530.1721(3 ) + projected semi - major axis of orbit , @xmath31 ( lt - s ) & 3.36671463(8 ) & 3.36671464(8 ) + longitude of periastron , @xmath64 ( deg ) & 1.35(2 ) & 1.36(2 ) + orbital eccentricity , @xmath65 & 1.91800(14)@xmath57 & 1.91796(15)@xmath57 + first derivative of orbital period , @xmath66 & 3.724(6)@xmath67 & 3.724(6)@xmath67 + first derivative of @xmath31 , @xmath68 ( @xmath69 ) & 1(2)@xmath54 & 1(2)@xmath54 + periastron advance , @xmath70 ( deg / yr ) & 0.0150(12 ) & 0.0150(13 ) + companion mass , @xmath71 ( @xmath72 ) & 0.223(14 ) & 0.223(15 ) + longitude of ascending node , @xmath73 ( degrees ) & 208.0(12 ) & 208.3(13 ) + orbital inclination angle , @xmath74 ( degrees ) & 137.1(8 ) & 137.3(8 ) + [ table:0437 ]     +   +    @xmath75{ngnoise.pdf } & \includegraphics[width=100mm]{gaussnoise.pdf }   \\ \hspace{-1.5 cm } \includegraphics[width=100mm]{redamp.pdf } & \hspace{-1.5 cm } \includegraphics[width=100mm]{redspec.pdf }   \\ \hspace{-1.5 cm } \includegraphics[width=100mm]{dmamp.pdf } & \hspace{-1.5 cm } \includegraphics[width=100mm]{dmspec.pdf }   \\ % \vspace{-1.5 cm } \end{array}$ ]    @xmath76{ngterms.pdf } \\ \includegraphics[width=120mm]{ngpdfs.png }   \\ % \end{array}$ ]    @xmath77{j043710cmnormres_scissored.pdf } & \includegraphics[width=85 mm , height=65mm]{10cmnormreshistogram.png } \\ \includegraphics[width=80mm]{j043720cmnormres_scissored.pdf }   & \includegraphics[width=85 mm , height=65mm]{20cmnormreshistogram.png }   \\ \includegraphics[width=80mm]{j043750cmnormres_scissored.pdf }   & \includegraphics[width=85 mm , height=65mm]{50cmnormreshistogram.png }   \\ % \end{array}$ ] in this paper we have introduced a method of performing a robust bayesian analysis of non - gaussianity present in the residuals in pulsar timing analysis , simultaneously with the pulsar timing model , and additional stochastic parameters such as those describing the red noise , and dispersion measure variations present in the data . deviations from gaussianity are described using a set of parameters @xmath37 that act to modify the probability density of the noise , such that @xmath78 describes gaussian noise , and any non zero values provide support for non - gaussian behaviour . the advantage of this method is that one may use a finite set of non - zero @xmath79 to model the non - gaussianity , without mathematical inconsistency . any truncation of the series still yields a proper distribution , in contrast to the more commonly used edgeworth expansion ( e.g. @xcite ) . we applied this method to two simulated datasets . in simulation one the noise was drawn from a non - gaussian distribution , and in simulation 2 it was purely gaussian . in simulation 1 , the effect of the non - gaussianity was to introduce a higher proportion of outliers relative to a gaussian distribution . this resulted in an overestimation of the toa uncertainties when assuming a gaussian likelihood , and decreased the precision with which the timing model parameters could be extracted compared to an analysis that correctly incorporated the non - gaussian behaviour on the noise .    in the second case we showed that the parameter estimates of the timing model parameters of interest were consistent when including , or not , the @xmath37 parameters , as is to be expected when the noise is gaussian . we then applied this method to the publicly available parkes pulsar timing array ( ppta ) data release 1 dataset for the binary pulsar j0437@xmath04715 . we detect a significant non - gaussian component in the non - thermal component of the uncorrelated noise , however as the non - gaussianity is most dominant in the lowest frequency data the impact on the timing precision in the pulsar is minimal , with only the parameter estimates of the power law dispersion measure variations being visible changed between the gaussian and non - gaussian analysis . janssen g. h. , stappers b. w. , kramer m. , purver m. , jessner a. , cognard i. , 2008 , in bassa c.,wang z. , cumming a. , kaspiv.m . , eds , aip conf . proc . 983 , 40 years of pulsars : millisecond pulsars , magnetars and more . phys . , new york , p. 633 skilling j. , 2004 , in fischer r. , preuss r. , von toussaint u. , eds , aip conf . 735 , bayesian inference and maximum entropy methods in science and engineering . , new york , p. 395
we introduce a method for performing a robust bayesian analysis of non - gaussianity present in pulsar timing data , simultaneously with the pulsar timing model , and additional stochastic parameters such as those describing red spin noise and dispersion measure variations . 
 the parameters used to define the presence of non - gaussianity are zero for gaussian processes , giving a simple method of defining the strength of non - gaussian behaviour . 
 we use simulations to show that assuming gaussian statistics when the noise in the data is drawn from a non - gaussian distribution can significantly increase the uncertainties associated with the pulsar timing model parameters . 
 we then apply the method to the publicly available 15 year parkes pulsar timing array data release 1 dataset for the binary pulsar j0437@xmath04715 . in this analysis we present a significant detection of non - gaussianity in the uncorrelated non - thermal noise , but we find that it does not yet impact the timing model or stochastic parameter estimates significantly compared to analysis performed assuming gaussian statistics . 
 the methods presented are , however , shown to be of immediate practical use for current european pulsar timing array ( epta ) and international pulsar timing array ( ipta ) datasets . 
 [ firstpage ]    methods : data analysis , pulsars : general , pulsars : individual