{'fraction_non_alphanumeric': 0.02900522655608829, 'fraction_numerical': 0.014729385339726146, 'mean_word_length': 4.569280731296608, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mass-Radius relationship for M dwarf exoplanets: Comparing nonparametric and parametric methods

M dwarfs, though the most abundant star in the galaxy, form only a small subset of stellar hosts with exoplanets with measured radii and masses. In this paper we analyze the Mass-Radius (M-R) relationship of planets around M dwarfs using M-R measurements for 24 exoplanets. In particular, we apply both parametric and nonparametric models and compare the two different fitting methods. We also use these methods to compare the results of the M dwarf M-R relationship with that from the Kepler sample. Using the nonparametric method, we find that the predicted masses for the smallest and largest planets around M dwarfs are smaller than that given by the relation fit to the Kepler data, but that the distribution of masses for 3 Earth Radii planets does not substantially differ between the two datasets. With future additions to the M dwarf M-R relation from TESS and instruments like the Habitable Zone Planet Finder, we will be able to characterize these differences in more detail, which will help illuminate the process of planetary formation and evolution around these stars. We release a publicly available Python code called MRExo which uses the nonparametric algorithm introduced by Ning et al. (2018) to fit the M-R relationship. Such a nonparametric fit does not assume an underlying power law fit to the measurements and hence can be used on samples spanning large mass and radii ranges. Therefore by not assuming a functional form, the fit is less biased. MRExo also offers a tool to predict mass from radius posteriors, and vice versa. This functionality can help inform observational strategies for radial velocity campaigns, such as TESS follow-up studies, as well as predict radii with microlensing planet masses.


INTRODUCTION
M dwarfs are the most common type of star in the galaxy (∼ 75%; Henry et al. 2006 Ricker et al. 2014) the hope is that we will soon discover hundreds of exoplanet candidates around them. While the discovery of these planets is interesting in itself, the comparison between them and the Kepler planets provide insight into the differing formation pathways of planets around M and FGK stars (Lissauer 2007). For example, the pre-main sequence lifetime of the star varies with its mass, which has an effect on the planetary migration process during its formation. Young M dwarfs are extremely active and exhibit high intensity XUV radiation which affects the inner planets and can potentially strip away the atmospheres of gaseous planets to leave rocky cores (Owen & Lai 2018;Owen & Wu 2017, show this for Kepler planets). There is also empirical evidence from Kepler (Dressing & Charbonneau 2015;Gaidos et al. 2016) and radial velocity (RV) surveys (Bonfils et al. 2013) which suggests that mass and radius distribution of planets is not identical for M and FGK dwarfs. These suggest that different physical processes may be at play, and pose a number of questions: Do smaller planets around M dwarf have more rocky compositions (Mulders et al. 2015a,b)? Is planet formation more efficient around M dwarfs (Dressing & Charbonneau 2015;Ballard & Johnson 2016;Ballard 2018)? If so, how does planet formation impact exoplanet chemical composition?
Probabilistic M-R relationships provide us with an empirical window into these questions, as they are closely related to distributions of exoplanet compositions. They also have very practical uses. For example, efficient planning of TESS follow-up radial velocity observations require an estimate of the planetary mass given a planetary radius. Future microlensing space missions like WFIRST (Green et al. 2012) will produce hundreds of exoplanets with the inverse problem of having a mass but not a radius. M-R relationships can be used to predict one quantity from the other. For this purpose, one needs a model for the M-R relationship that best balances the trade-off between the prediction's variance (i.e. the width of the range of possible masses for a given planet) and its bias in the predicted values (i.e. the difference between the true mass and the mean predicted mass). M dwarfs will also offer potential targets for atmospheric characterization of Earth-like habitable zone planets using transit spectroscopy due to the lower stellar brightness and relatively large planet to star radius ratio. It has also been shown that the spot induced radial velocity jitter is reduced in near infrared (NIR) (Marchwinski et al. 2015), which is where the spectral energy distribution for these stars peak. For these reasons, many more M dwarf planets will have their masses measured in the near future by instruments such as HPF (Mahadevan et al. 2012), CARMENES (Quirrenbach et al. 2016), NIRPS (Wildi et al. 2017), IRD (Kotani et al. 2018), SPIRou (Artigau et al. 2014), iSHELL (Cale et al. 2018), GIANO (Claudi et al. 2017), and NEID (Schwab et al. 2016). Here we set the stage for these future datasets by assessing the dependence of M-R model choices on mass and radius predictions, especially as a function of stellar type. Substantial past efforts have been performed in studying M-R relations in recent years. A summary of this is presented in Ning et al. (2018). Several of the widely used M-R relationships include Weiss & Marcy (2014); Wolfgang et al. (2016); Chen & Kipping (2017), the latter of which also introduced a publicly available Python package called Forecaster. The M-R model underlying Forecaster uses a broken power law to fit the M-R relationship across a vast range of masses and radii, in recognition of potential changes in the physical mechanisms responsible for the planetary formation at different mass regimes. However, as has been discussed in Ning et al. (2018) such a restrictive parametric model can portray an incom-plete picture since we do not know the true functional form of the underlying relationship, such as whether it is a power law to begin with. Conversely, a nonparametric model offers more flexibility in the fit, which can be advantageous when the goal is to obtain predictions that best reflect the existing dataset. Ning et al. (2018) introduced a nonparametric model for the M-R relationship which uses Bernstein polynomials, a series of unnormalized Beta probability distributions. We apply this model to M dwarf exoplanets in preparation for future, larger datasets of small planets around small stars, which are less likely to be fully described by highly parametric models.
While adapting the methodology of Ning et al. (2018) to a new dataset, we have developed and are now offering a publicly available Python package called MREXo, inspired by the useful community tool Forecaster. MREXo can not only be used as a predicting and plotting tool for our M dwarf and Kepler M-R relationships, but can also be used to fit a M-R relationship to any other dataset (but see §6 for a discussion about the minimum dataset size at which a nonparametric fit becomes useful). This makes it a powerful tool for exoplanet population studies and a probe for potential differences in composition across samples. This package uses the open source tools of Python deployed with parallel processing for efficient computation. It also offers a fast predictive tool for the M dwarf and Kepler sample M-R fitsused in this paper, so that either mass or radius can be predicted from other.
The rest of the paper is structured as follows. In §2 we discuss the input dataset, and in §3 we discuss the parametric and nonparametric fitting process and algorithms followed. In §4, we describe the Python package MRExo that we release along with this paper. Then, in §5 we discuss the results of the fits, where we compare the parametric fit with the nonparametric and the M dwarf M-R relationship with that from the Kepler exoplanet sample. In §6, we explain the simulation performed to test the efficacy of the nonparametric method. We end with the discussion in §7 and conclude in §8.

DATASET
Fitting a M-R relationship requires a sample set with confirmed mass and radius measurements 1 The mass and radius values for exoplanets used in this work are obtained from the NASA Exoplanet Archive, which was last accessed on 14th December 2018 (Akeson et al. 2013). Figure 1 shows the 24 M-R data points that we have used, colour-coded by host star temperature. The mean values for mass and radius, along with their respective measurement uncertainties, are shown in Table A1 in the Appendix. We include mass values estimates from both RV and Transit Timing Variations (TTVs). This sample is hereafter referred to as the M dwarf dataset in this paper. We specifically chose to omit the three planets discovered by Direct Imaging with both radius and mass constraints, as these planets have substantially larger orbital separations than the other planets in our sample. Furthermore, the directly imaged planets have their masses and radii modelled and not directly measured.
To limit ourselves to the M dwarfs, we restrict our host star sample to T eff < 4000 K. To exclude brown dwarf companions, we restrict ourselves to planetary masses (M p < 10 M J ). We chose to use the most-recent values from the Exoplanet Archive rather than the default values. We manually verified the reported values from the Exoplanet Archive with the literature refer-1 The nonparametric framework employed here can not handle measurement upper limits, since the maximum likelihood estimation method used in Ning et al. (2018) does not allow censored data; hence planets with only upper limits are excluded from the results presented in §5.2. Adapting this methodology to include upper limits is an area for future work. ences. For the TRAPPIST-1 system we use the Grimm et al. (2018) values for the mass and radius measurements. Note that of the 24 exoplanet measurements in our sample, 7 belong to the TRAPPIST-1 system and thus our M-R joint distribution is heavily influenced by the TRAPPIST-1 system in the small planet (R < 2 R ⊕ ) regime. For the purposes of comparison we also use the Kepler dataset from Ning et al. (2018) to compare our M dwarf results with that from a larger Kepler sample, which consists of 127 M-R measurements and does not have any T eff restrictions on the host stars. This sample is hereafter referred to as the Kepler dataset in this paper.

FITTING THE DATASET
To fit these datasets we use two recently proposed methods: a parametric method proposed by Wolfgang et al. (2016), and a nonparametric method by Ning et al. (2018). The two methods and their M-R fits are detailed and compared in subsequent sections. We subsequently perform the following analysis: 1. Comparing the parametric and nonparametric results for the M dwarf sample.
2. Comparing the M dwarf and Kepler sample M-R relations using the nonparametric model.
3. Performing a simulation study to show the efficacy of the nonparametric method as a) the sample size increases, and b) the intrinsic astrophysical dispersion in the sample increases (assuming a power law with an intrinsic dispersion).
We compare the parametric method with the nonparametric technique because the planetary masses and radii of the M dwarf sample currently appear to resemble a power law. This may not remain the case as we accumulate more data for M dwarf planets; indeed, with over 100 mass and radius measurements, the Kepler dataset is not currently well described by a power-law model (Ning et al. 2018). Additionally, nonparametric methods offer less biased predictions, as long as the dataset it was fit to is representative of the underlying distribution. In contrast, nonparametric methods often perform poorly for small data sets and can easily overfit the data, while parametric methods are easier to implement and are valuable for gaining an intuitive understanding of the dataset. Acknowledging both the simplicity and easy insight from a parametric model and the improved predictive capabilities of a nonparametric model for larger datasets, we compute both and compare the outputs for the current M dwarf sample (see §5). In addition, we run simulations to assess at what dataset size does the nonparametric technique work well (we also test the fitting as the intrinsic dispersion in the dataset increases).
We note that the nonparametric method can not handle upper limits in its current implementation. In the TESS era of planet discoveries, we will soon have many exoplanet candidates with just upper limits, and so this is a clear area for future statistical development of this nonparametric approach. For the time being, we incorporate the upper limits in the parametric approach, which can accommodate censored data like upper limits.

Parametric method: Hierarchical Bayesian Modelling
The parametric method we use is heavily borrowed from the probabilistic hierarchical Bayesian model introduced in Wolfgang et al. (2016). This model consists of a power law to describe the mean planet mass as a function of radius, plus a normal distribution around that line to describe the intrinsic, astrophysical scatter in planetary masses at a specific radius. The model we employ here is very similar to Wolfgang et al. (2016, Equation 2), except that we use a lognormal distribution to describe the intrinsic scatter. This change causes the intrinsic dispersion to be constant in log space, which Chen & Kipping (2017) show is a better descriptor of the scatter over a larger mass range than was considered in Wolfgang et al. (2016). The parametric model used here is therefore: (1) where C is the normalization constant of the mean power law (once it is converted to the linear mass scale), the γ is the power-law index, and the σ M is the intrinsic scatter, in terms of log(M ). The ∼ symbol implies that the masses are drawn from a probabilistic distribution (as opposed to a '=' sign for a deterministic model); this distribution is what models the astrophysical scatter in planetary masses. Here it is parameterized as a Gaussian centered on a line, which when converted to linear scale is a log-normal distribution centered on a power law. Evaluating this model within a hierarchical framework allows the upper limits in our dataset to be incorporated into the inference of C, γ, and σ M , along with arbitrary measurement uncertainty for individual data.

The Nonparametric method -Bernstein polynomials
As explained in the introduction, we also employ the nonparametric model introduced by Ning et al. (2018) to fit M-R relationships, with an eye toward future M dwarf dataset sizes which will likely benefit from a more flexible approach. The nonparametric method fits the joint M-R distribution, and hence can treat either mass or radius as the independent variable, and can be used to predict one from the other. This method uses a sequence of Bernstein polynomials as the basis functions to fit a nonparametric M-R relationship; when normalized these Bernstein polynomials are identical to Beta probability distributions. Hence, the model is equivalent to a mixture of unnormalized Beta probability distributions. We translate the code presented in Ning et al. (2018) from R to Python and release it in a publicly available package called MRExo; this code is further discussed in §4.
The nonparametric approach uses these Bernstein polynomials as the basis functions to fit the joint distribution of masses and radii f(m,r). This is detailed in §2 of Ning et al. (2018): where w kl is the kl-th element of the matrix w, which is a set of weights corresponding to individual Bernstein polynomials. Each weight is positive, and the sum of individual weights equals unity. The weights describe how much each term in the series contributes to the joint distribution. The degrees for the polynomials in the mass and radius dimensions are represented by d and d , respectively. Here m and r depict the continuous variables which represent mass and radii. M , M , R and R represent the lower and upper bounds in mass and radius, respectively. To fit for the weights (w kl ), we use maximum likelihood estimation (MLE). The likelihood (Equation 4) includes measurement uncertainties as normal distributions, hence Equation 4 modifies Equation 2 to introduce a convolution of the normal and Bernstein polynomials.
The measured values for mass and radii are assumed to be drawn from a normal distribution centered on the true value, with a standard deviation equal to the measurement uncertainty. This produces the joint distribution shown in Equation 4, where M obs i , σ obs M i are the mass observations and their uncertainties and likewise R obs i , σ obs R i for radius. N denotes a normal distribution. After optimization for the weights using the MLE method via the Python package SciPy (Oliphant 2007), the joint distribution for the M dwarf dataset can be plotted (see Figure 3

b).
A key consideration in fitting this nonparametric model to a dataset is identifying the optimum degree for the Bernstein polynomial series. We use the cross validation method as explained in Ning et al. (2018) to find that the optimum values for d and d are both 17. In total there are d× d weights which serve as our 'parameters' to fit for. While 17 2 = 289 parameters may seem excessive for a dataset size of 24 planets, in practice there are only five non-zero weights in the series (see Figure 2 for a pictoral representation of this). This highlights one of the advantages of using Bernstein polynomials as our basis functions: estimating the weights in this series is self-regularizing, meaning that the smallest number of non-zero coefficients is automatically found. Additionally, using Bernstein polynomials efficiently reduces the number of nonzero free parameters; if one instead used a mixture of Gaussians to fit the joint distribution, that fit would require at least 3 times as many free parameters (amplitude, mean, and standard deviation, per component, rather than just amplitude).
After finding the weights via MLE, we repeat the process using the bootstrap method, which estimates the uncertainty in this fit. We do this by resampling the dataset with replacement and running the fitting routine again for each bootstrap. In regions without data points, the Bernstein polynomials revert to the overall mean of the distribution.
We found that the Bernstein polynomials may behave counter-intuitively at the boundaries of the joint distribution (see Appendix A in Ning et al. (2018)). To address this issue, we fix the first and last row and column of Bernstein polynomials to have zero weights (w 0j , w dj , w i0 , w id ), and as such they do not contribute to the joint distribution. This was not done by Ning et al. f (m, r|w , d, d , M obs (2018) since their Kepler sample had enough samples near the boundaries to constrain the fit.

MREXO
In this section, we shall discuss an important contribution made in this paper. We translate Ning et al. (2018)'s R code 2 into a publicly available Python package called MRExo 3 . This can be easily installed using PyPI and has extensive documentation and tutorials to make it easy to use.
MRExo offers tools for fitting the M-R relationship to a given data set. In this package we use a cross validation technique to optimize for the number of degrees. We then fit the joint distribution ( § 3.2) to the sample set; this can then be marginalized to obtain the conditional distribution, which we can use to predict one variable from the other. We bootstrap our fitting procedure to estimate the uncertainties in the mean M-R relation. Further, MRExo is equipped with dedicated and easy to use functions to plot the best-fit conditional M-R relationships as well as the joint M-R distribution 4 . Crucially, MRExo also predicts mass from radius, and radius from mass. For example, in the case of planets discovered using the transit method, the feasibility of an RV follow-up campaign can be evaluated by predicting the estimated mass and its confidence intervals given the measured radius and its uncertainty. Another feature of this package is that it can accommodate radius (or mass) posterior samples from separate analyses, which are then used to compute the posterior predictive distribution for mass . Along with the MRExo installation, the results from the M dwarf sample dataset from this paper and the Kepler exoplanet sample from (Ning et al. 2018) are included.
The degree of the Bernstein polynomials (d) approximately scales with the sample size (N ). Since the number of weights goes as d 2 , the computation time involved in the fitting a new M-R can soon start to become prohibitive. Therefore we also parallelize the fitting procedure and the bootstrapping algorithm. As an example, the M dwarf sample (N = 24; d = 17) took about 2 minutes to perform cross validation, fit a relationship, and do 100 bootstraps on a cluster node with 24 cores and 2.2 GHz processors. The Kepler sample (N = 127; d = 55) took about 36 hours for the cross validation, fitting, and 48 bootstraps. We realize that the fitting computation time would start to become prohibitive as the sample size increases 200; therefore we plan to optimize the code further by benchmarking, optimizing floating point operations, and correcting the precision requirements in the integration step. However, this time intensive step of cross validation and fitting is only nec-essary if the user needs to fit their own relationship. To run the prediction routine on the pre-existing M dwarf or Kepler samples is fairly quick and takes less than a second for a prediction. In order to do a large number of predictions as part of a larger pipeline or simulation, the user can also generate a look-up table which makes the calculations even faster (the function to generate and use the look-up table is provided with the package).

Parametric fit results
As discussed in §3, there is a trade-off between the flexibility and lower bias of nonparametric models and the interpretability and lower predictive variance of parametric models. At present, the M dwarf dataset only consists of 24 planets, and their masses and radii seem to be well approximated by a power law. As such, we fit the parametric model described by Equation 1 to the same M dwarf dataset, to serve as a basis for comparison to the nonparametric results. Currently the fit looks to be fairly reasonable by eye (see Figure 3e), but this is not guaranteed to be the case as more measurements are obtained (indeed, it is clear from the fit to the Kepler data, which consist of > 100 planets, that a single power law is not a sufficient model for that data; see Figure 3f).
The best-fit parameter values for the model described in Equation 1 are C = −0.130 +0.081 −0.055 , γ = 2.14 +0.11 −0.16 , and σ M = 0.184 +0.077 −0.021 ; these values were found by identifying the mode of the joint 3-dimensional hyperparameter posterior, and the error bars represent the marginal central 68% credible intervals. This power-law slope is steeper than that found by Wolfgang et al. (2016) for nearly all of the datasets they consider (most lie within 1.3 < γ < 1.8). There are two possible explanations for this. First, the dataset used by Wolfgang et al. (2016) and, in fact, nearly all previous M-R results except Neil & Rogers (2018) -is dominated by FGK dwarf planet hosts. This result could therefore be driven by intrinsic differences between the planets that form around Sun-like stars and those that form around M dwarfs. However, Wolfgang et al. (2016) also use only a limited radius range for their dataset (R < 4 R ⊕ or R < 8 R ⊕ ), as their focus was on super-Earths. As shown by Ning et al. (2018), the slope for 1 < R < 4 R ⊕ is shallower than that for 4 < R < 10 R ⊕ , and so one would expect that a single power-law fit to the entire radius range would result in a steeper slope than that fit just to R < 4 R ⊕ . This is corroborated by the fact that the slopes fit to the R < 8 R ⊕ radius range are steeper than those fit to the R < 4 R ⊕ range. Additionally, Chen & Kipping (2017) find that the segment of their broken power-law relation which most closely matches radius range of our M dwarf dataset (1 R 11 R ⊕ ) has a power-law slope of 1 0.59 = 1.7. We note that the fact that there is a noticeable difference between the slopes fit to different radius ranges is evidence that a more flexible relation such as the non-parametric one developed in Ning et al. (2018) is needed. Another notable difference between these two relations is the predicted mass for a 1 R ⊕ planet. This information is contained in C, the mean power-law normalization constant. In linear units, the best fit C for the M dwarf dataset is 10 −0.13 = 0.74 M ⊕ . This is significantly lower than the power-law constants found in Wolfgang et al. (2016): it is inconsistent with the R < 8 R ⊕ fit (C = 1.5 M ⊕ ) at 4σ, and with the R < 4 R ⊕ fit (C = 2.7 M ⊕ ) at 7σ. At least two effects are contributing to this difference. First, very few planets around FGK dwarfs with R < 1 R ⊕ have had their masses measured and constrained to be M < 1.5 M ⊕ . Because of the publication bias quantified in Burt et al. (2018), wherein planetary masses are required to reach a certain significance threshold to be published, small, low-mass planets are more likely to be left out of the FGK dwarf datasets. This causes the relation to be fit to the more massive small planets that do appear in the literature, which in turn causes C to be high. Conversely, the M dwarf dataset is strongly affected by the presence of the TRAPPIST-1 planets (Gillon et al. 2017), whose masses are reported by Grimm et al. (2018) to both small and quite precise. This one analysis of a single system dominates the M dwarf dataset at 0.7 < R < 1.2 R ⊕ and therefore the best-fit value of C -a caveat that all potential users of these M dwarf relations should keep in mind.
To determine whether the differences highlighted above are due to the differences in the host star population or due to differences in the models or the considered radius ranges, we also fit this simple parametric model to the Kepler dataset from Ning et al. (2018). While we don't expect this model to be a good fit to the data, it offers an apples-to-apples comparison to the above M dwarf parametric results, and its shortcomings highlight the advantages of more flexible nonparametric models. The best-fit parameter values for the Kepler dataset are C = −0.0250 +0.093 −0.111 , γ = 2.13 +0.11 −0.14 , and σ M = 0.457 +0.026 −0.044 . Therefore there is no statistically significant difference between the M dwarf planetary mean power-law and the one that would be fit to the planets around FGK dwarfs. More specifically, the two slopes are nearly identical; while the Kepler normalization constant is larger than the M dwarf constant, the difference is at the 1σ level. This comparison is consistent with the result of Neil & Rogers (2018), which used a sample size of 6 M dwarf planets. On the other hand, there is a statistically significant difference in the intrinsic scatter for the M dwarf and Kepler datasets, with the Kepler dataset having more variation in planet mass at a given radius. However, this may be at least partially affected by the differ-ences in the sample sizes (N = 24 vs. N = 127) and the fact that it takes more time to get significant mass measurements for planets at the low mass side of the intrinsic scatter distribution. As we obtain more masses for transiting M dwarf planets with instruments like HPF, we will be able to test whether the difference in the intrinsic scatter remains statistically significant.
All told, the comparison between the M dwarf and Kepler parametric model fits does not yield very much insight into the differences between the planetary populations. Because the parametric model is not able to capture detailed features of the M-R relation, it may hide some important differences between the two datasets once the M dwarf dataset becomes large enough to warrant detailed analyses of these features. Looking ahead to the mass-radius dataset we will have by the end of the TESS mission, we apply the nonparametric model to the current M dwarf dataset and perform the comparison between it and the nonparametric fits to the Kepler dataset to serve as a basis for comparison to this future work (see §5.4).

Nonparametric fit
We fit the M dwarf exoplanet dataset shown in Figure 1 using the nonparametric approach described in §3.2; the results of this fit is displayed in Figure 3 a,b. Using the cross validation method the optimum number of degrees selected are both equal to 17, giving 17 2 = 289 total weights. We note that our algorithm automatically forces most of the weights to zero to prevent over-fitting, especially for small datasets. For this dataset, only five weights are non-zero (see Figure 2).
We also plot the conditional distributions f (m|r) and f (r|m), where the conditional distribution f (m|r) is the ratio of the joint distribution (Equation 2) to the marginal distribution of the radius f (r); likewise for f (r|m)(see Equations 1 and 2 of Ning et al. (2018)). The conditional distribution spans regions of M-R Figure 2. The optimal weights chosen by cross validation for the M dwarf sample. The cross validation procedure selected 17 degrees as the optimum for both mass and radius which gives 17 2 = 289 weights, of which only 5 (4 are visible, the 5th one at 9,12 is difficult to see in this stretch) are non zero.
space where there is no data. This is due to the symmetrical nature of the Bernstein polynomials and indicates that the current dataset is not yet large enough to be effectively described by a nonparametric method; this is further discussed in §6.
From Figure 3b we note that the red conditional distribution-f (m|r) is not the same as the blue conditional distribution-f (r|m). This is because when the joint distribution (see Figure 3 b,d) is marginalized to obtain the conditional, the distribution behaves differently for each axis. The two conditional distributions would be the same for a dataset with very little error and symmetric, equal scatter along the entire M-R relation (see §6 for an illustration of this). The Kepler dataset effectively illustrates how localized areas of large scatter in one dimension can drive differences between f (m|r) and f (r|m). In particular, at log(m) ∼ 0.5 there is large scatter in radius, which is fit by the mean M-R relation in f (m|r) but is rep-resented by the distribution (the width of the light shaded region) around that mean relation in f (r|m). While both conditionals were computed from the same joint distribution (Figure 3d), asymmetries such as this in the joint distribution can produce conditionals that qualitatively look very different. We also note that truncating the probability distributions for integration (that is, using finite bounds to compute f (r) and f (m) for the conditionals) also contributes to some differences between the two conditional distributions. This can be seen further in the simulation in Figure 5, where the disparity between the two conditional distributions increases as the dispersion from the power law increases. Therefore, to predict one quantity from the other it is imperative to use the right conditional distribution.
A common and fair concern about nonparametric methods is their ready ability to overfit data. The cross-validation method we adopt to choose d and d is designed to minimize overfitting while maximizing the predictive accuracy of the M-R relation. Both underfitting and overfitting produce high predictive error (that is, the predicted value is far from the true value); because the cross-validation method identifies the degree which minimizes the predictive error, it finds the optimum d that balances the trade-off between the two. That said, we acknowledge that by eye there are wiggles in the M dwarf M-R relation that don't appear to be supported by the dataset. This is likely due to the sparsity of the dataset in certain radius and mass ranges; where there is no data, the relation tends toward the mean of the closest nonzero term. This is why we perform the simulation study in §6, to find at what size dataset does our non-parametric model effectively describe the conditionals. We plot the conditional probability density functions (CPDFs) of planetary masses for planets with R = 1, 3, and 10 R ⊕ in Figure 4. These CPDFs show the distribution of planetary masses that would be predicted for planets at those radii. Assessing the differences visually, the parametric and nonparametric M dwarf fits have the same median, but different spreads. This difference is expected due to the bias vs. variance trade-off in parametric vs. nonparametric modeling: if the parametric model is a poor description of the data, the nonparametric fit will be less biased but have a higher variance than the parametric estimator. Because the medians of the CPDFs are similar between the non-parametric and parametric models, it appears that, with this current dataset, a powerlaw is a decent fit to the data. That said, we emphasize that this will likely not remain the case as the M dwarf dataset grows, as we have seen for the Kepler dataset.

Comparing Nonparametric vs. Parametric fitting
To quantify the differences between the mass predictions produced by these two model choices, we perform the two-sample Kolmogorov Smirnov (KS;Kolmogorov 1933;Smirnov 1948) and Anderson Darling (AD; Anderson & Darling 1952;Scholz & Stephens 1987) tests on samples drawn from these CPDFs (see Table 1). These statistical tests assess whether we can reject the null hypothesis that two datasets were drawn from the same distribution, with the KS test being more sensitive to the location of the median and the AD tests more sensitive to differences of the tails.
Importantly, the reported significance of these tests depends on the size of the datasets, with small differences between the samples becoming more significant as the datasets grow. Since we have the functional form of these CPDFs from our fitting procedure, we have a choice in how many samples we can draw from them, and therefore how significant we report them to be 5 . As such, the results in Tables 1 and 2 should be interpreted on a comparative basis, not on an absolute basis: rather than concluding that the CPDFs for 1 R ⊕ planets are different on a 5 For every comparison in Tables 1 and 2, we use 100 samples from each CPDF. We also note that the assignment of a p-value to a given statistic depends on a theoretical statistical distribution for the variation in that statistic given randomly drawn datasets from the same distribution. To test this theoretically assigned pvalue, we generated reference distributions of KS and AD statistics based on 20,000 randomly drawn datasets (each of which also contain N = 100 planets) from the CPDFs, compared against themselves for a total of 10,000 comparisons (and 10,000 values for the KS and AD statistics). We then identified at what quantile is the median of the nonparametric vs. parametric KS statistic distribution with respect to the distribution of KS statistics produced by comparing the M dwarf nonparametric CPDFs against themselves. These values are what is reported in Table 1 (and an analogous comparison but between the Kepler and M dwarf nonparametric CPDFs with respect to the self-compared nonparametric M dwarf CPDFs is reported in Table 2). It turns out that this more careful approach is in good agreement with the theoretical p-values provided by the R package kSamples, built from the work of Scholz & Stephens (1987). If the median of the nonparametric vs. parametric KS statistic distribution fell completely outside that of the nonparametric self-compared KS distribution, the p-value is reported as < 10 −4 , as we only generate 10,000 pairs of datasets for these KS and AD statistic distributions. Figure 4. The conditional probability density functions CPDFs of mass given radius for planets at 1, 3, and 10 R ⊕ ; these show the predicted mass distributions at these radii from the nonparametric fit to the M dwarf sample (top); the nonparametric fit to the Kepler sample (middle); and the parametric fit to the M dwarf sample (bottom). The means of both M dwarf fits are at lower mass for 1 and 10 R ⊕ planets than for the Kepler fit, indicating that the smallest and largest planets are less massive around M dwarfs than around Sun-like stars. Interestingly, the means of the predictive masses for 3 R ⊕ planets are similar. . (Colour version of this figure is available in the online journal.) statistically significant level while the CPDFs for 10 R ⊕ are not, the correct interpretation is that it will take a smaller number of mass measurements of planets at 1 R ⊕ to distinguish between the parametric and nonparametric models than it will for planets at 10 R ⊕ (see Table 1). Along the same vein, an even smaller number of mass measurements would be needed to distinguish between the parametric and nonparametric models for 3 R ⊕ planets. Because we only have 24 planets in our M dwarf dataset in total, let alone at any single radius, these model comparisons will need to be reassessed in the future with more data.

Comparing M dwarf M-R vs. Kepler M-R
An important point of this paper is to compare the M-R relationship from a Kepler exoplanet sample to an M dwarf exoplanet sample. As discussed in §1, there is empirical evi-dence that the planetary radius and mass distributions differ between M and FGK dwarfs. This is further illustrated empirically by our fit to the M dwarf and Kepler sample (Figure 3): there are visual differences between the joint and conditional distributions for the M dwarfs and Kepler. Additionally, Figure 4 shows that the median value of the CPDFs of Kepler and M dwarf masses do not coincide for the smallest and largest planets; this is true whether the nonparametric M dwarf relation or the parametric M dwarf relation is used. Therefore, if the Kepler M-R relationship was used to predict the mass of a transiting M dwarf exoplanet like TRAPPIST-1 b, it would produce mass predictions that were too large on average. Conversely if used for a non transiting planet like Proxima-b to predict its radius from its mass, the prediction would be erroneous. This illustrates that the conditional density functions are different for the two samples. We also perform the KS and AD tests to quantify differences between the M dwarf and Kepler CPDFs (see §5.3 for discussion about these tests). We find that it will take << 100 planets to distinguish between the M dwarf and Kepler mass predictions at 1 and 10 R ⊕ , while even a dataset of 100 3 R ⊕ planets does not illuminate any statistically significant difference between the mass predictions. Therefore, the M dwarf M-R relation could be most different from the FGK-dwarf M-R relation at the smallest and the largest planet sizes.

SIMULATION STUDY USING MREXO
As discussed in §3, we run a simulation to test the effectiveness of the nonparametric method as a fitting and predictive tool. In particular, we visually assess the ability of the Bernstein polynomial model to qualitatively reproduce a known underlying distribution as a function of sample size and astrophysical scatter. We set the known distribution to be a power law in M-R space with a slope of ∼ 2.3 (mass as a function of radius). We simulated synthetic datasets We note that in the simulation the nonparametric algorithm fits the simulated dataset well in the case of low intrinsic dispersion even for small datasets (∼ 20 points; see Figure 6). However, as we increase the intrinsic dispersion, the fit compensates by increased uncertainties and visual departures from the known underlying distribution see Figure 5). Therefore, this simulation confirms that the utility of the nonparametric method increases as the sample set increases in size ( Figure 6). Importantly, we note that this simulation uses a simple power law, for which a parametric fit would be a sufficient descriptor of the data. However, we do not know the true underlying distribution of nature's exoplanet masses and radii, and so we are justified in applying a nonparametric model to real data with a large enough sample size.  sample. A parametric model is useful for small datasets where it offers an easy means to develop an intuitive understanding for the data. It can also be used when there is a known physical process that is being fit or tested for. However, in cases where we need to explore the relationship and develop a forecasting or predicting routine, nonparametric methods offer a more flexible and unbiased option. Furthermore, in cases where the underlying distributions of the variables is unknown, or the data could have small scale structure which can not be captured by a power law, the nonparametric method is a better predictor. For this M dwarf sample, the parametric and nonparametric methods offer similar median predictions in the realm of limited sample size (see Figure 4). However, as has been shown by Ning et al. (2018) with the Kepler sample, a simple parametric power law (or broken power law) does not necessarily Figure 6. The nonparametric fit to simulated samples of 20, 50, and 100 planets with the same amount of intrinsic scatter. The fit improves as the sample size increases. One can also note that the conditional f(m|r) and f(r|m) distributions begin to converge as the sample size increases. (Colour version of this figure is available in the online journal.) capture the features in the data as the sample size grows. A parametric predicting method will not necessarily reflect the data with the same accuracy and has the potential for higher bias, meaning that the predicted masses or radii will be farther from the true masses or radii. Therefore we propose this nonparametric method for the M-R relationship and offer the Python pack-age MRExo as a practical tool for making these predictions.

M dwarf vs Kepler
Neil & Rogers (2018) fit a M-R relationship on a sample of M dwarf exoplanets and compare it to a sample of exoplanets orbiting FGK dwarfs, to try to find evidence for any dependence on the host star mass. They do so by expanding upon the HBM approach introduced in Wolfgang et al. (2016), to include another dimension that accounts for a possible host star mass dependence. As their M dwarf sample consisted of only 6 planets with radii and mass measurements, they find that their M-R power law was slightly shallower than that for their FGK sample; yet still consistent with no host star mass dependence. Our sample size of 24 planets agrees with this when fit with a similar parametric relationship.
On the other hand, we find in §5.4 that the Kepler and M dwarf predictive mass distributions are different for the smallest and largest planets, in the sense that the smallest and largest planets are on average less massive around M dwarfs than around Sun-like stars. This result holds for both the nonparametric and the parametric M dwarf fit. The difference appears in these comparisons but not in the parametric comparison. This is because as shown by Ning et al. (2018) the parametric model is not flexible enough to fit the larger Kepler dataset, and hence does not provide a reasonable baseline for comparison.
Astrophysically, this could indicate that the low-mass protoplanetary disks around M dwarfs result in lower-mass gas giants, which would explain the lower predictive masses for the 10 R ⊕ . For the smaller planets that are likely terrestrial, a lower mass at a given radius means a higher silicon-to-iron ratio, which may suggest that M dwarf disks process their refractory elements differently than Sun-like protoplanetary disks. Interestingly, the predictive mass distributions for 3 R ⊕ planets are similar between M-dwarf and FGK-dwarf hosts. This could indicate that the process of sub-Neptune formation is independent of the mass of the protoplanetary disk. That said, the predictive mass distributions between the M dwarf and FGK planetary samples do overlap substantially, and the M dwarf dataset is still quite small and is dominated by the TRAPPIST-1 planets. There-fore, these preliminary results should be revisited once more transiting M-dwarf planets are discovered and followed up with ground-based observations.

Future Prospects
As was seen with the Kepler sample, as the sample size increases we can start to unearth interesting phenomena that would be hard to discern from small data sets. In addition, more precise measurements with smaller error bars would help in more accurate M-R model fits. This would further unveil features and regions that were indistinguishable earlier. This will particularly hold true as the M dwarf M-R space starts to fill up with radii from TESS (Ricker et al. 2014

CONCLUSION
In this paper we fit the mass-radius relationship for a sample of 24 exoplanets around M dwarfs using nonparametric and parametric methods. Considering the small sample size, the measurements are currently well described by a power law, which we fit with a parametric hierarchical Bayesian model. We then discuss the deficiencies of parametric methods, and the utility of nonparametric models, especially as the sample size increases. To further illustrate this point we also run a simulation study showing how the nonparametric fit changes with a change in sample size and dispersion in the sample. We then compare the nonparametric and parametric results, finding them to be similar on average but currently with a larger variance for the nonparametric mass predictions. We also discuss differences in the M-R relationship for M dwarf vs Kepler exoplanets and note that the predicted conditional probability density functions differ for the smallest and largest planets. Using this comparison we illustrate the utility of an M dwarf M-R relationship in era of the exciting new discoveries with TESS and groundbased precision RV instrumentation.
We also introduce a new Python package called MRExo which can be used as a predictive tool, as well as to fit nonparametric M-R relationships to new datasets. This code is available to the community to use in their own applications.