An Automatic Deconvolution Method for Modified Gaussian Model using the Exchange Monte Carlo Method: Application to Reflectance Spectra of Synthetic Clinopyroxene

Copyright: © 2016 Hong PK, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Introduction
Remote sensing of reflectance spectra of Earth and other planetary bodies can be useful for identifying mineral distribution on their surfaces, especially in remote regions that are exceedingly challenging to perform field-based investigation, and those planetary surfaces yet to have in situ observation, mapping, characterization, sampling, and analyses [1][2][3][4][5][6][7][8]. Many different factors, however, can influence the surface spectra, such as various alteration and weathering processes, and observational conditions [2,3,5,[9][10][11][12]. Because a reflectance spectrum is a complex non-linear mixture of the above mentioned factors [13][14][15][16], it is highly challenging to segregate each factor and extract the true mineral spectra, based solely on remotely observed reflectance spectra, and thus confidence in the resulting signatures should be gained by comparing with the reference spectra obtained by field or laboratory measurements [17]. Though challenging, investigating the spectral change due to the variation of elemental composition should not be avoided, since it is one of the ultimate goals of remote sensing of reflectance spectra of planetary surfaces [2,18]. Compared to terrestrial surfaces, those of extraterrestrial bodies such as the Moon and asteroids are not covered by liquid water and vegetation, and have negligible to no atmosphere, and thus may be considered to be the best places to observe the true nature of mineral spectra [19][20][21][22][23][24][25]. Yet, there are many factors to contaminate reflectance spectra of such planetary surfaces such as regolith particles, space weathering, and horizontal and vertical mixing by impact cratering, thus identifying the variation of mineral distribution on extraterrestrial bodies is still difficult [10,22,[26][27][28][29][30]. Therefore, studying the spectral change using simple pure minerals that compose Earth and planetary surfaces is a critical foundational step to analyze reflectance spectra and is a prerequisite of reflectance spectroscopy, in light of continued application to planetary surfaces [8,[31][32][33][34][35].
Deconvolution of reflectance spectra has been a common procedure for interpreting the experimental data of minerals and observation data of terrestrial and extraterrestrial surfaces. Among the most recognized spectral deconvolution methods in planetary science is the modified Gaussian model (MGM) [36,37]. In their paper, Sunshine et al. [36] showed that for 1 µm absorption band of orthopyroxene, Gaussian functions in the wavelength space fit better than Gaussians in the frequency space. The MGM express reflectance spectrum by the following function: where R is reflectance, x is wavelength, and s k , µ k and σ k are the strength (amplitude), center (mean) and width (standard deviation)

Abstract
Deconvolution analysis of reflectance spectra has been a useful method to infer mineral composition and crystal structure. Many of the recent deconvolution analyses of reflectance spectra of major rock-forming minerals, such as olivine and pyroxene, have been based on a modified Gaussian Model (MGM). The numerical algorithm of the widely used MGM, however, utilizes the steepest descent method, which has a local minima problem. With inaccurate initial parameters, the steepest descent method converges into a local minimum, thus the analyzer must manually adjust initial parameters and calculate the model repeatedly to obtain the desired solution. In order to avoid the local minimum problem, we utilized Bayesian spectral deconvolution with the exchange Monte Carlo method, which is an improved algorithm of the Markov chain Monte Carlo method, aimed to both avoid local minima traps and remove the arbitrariness originated from initial parameters. We applied the model to visible to near infrared reflectance spectra of 31 synthetic clinopyroxene samples with wide ranging Mg, Fe and Ca compositions (solid solution). We obtained results consistent with the previous studies based on conventional MGM analyses, suggesting that the exchange Monte Carlo method can yield results consistent with the conventional MGM analyses purely based on the observed data. We also find that the center wavelengths of 1 µm absorption bands of high-Ca pyroxene samples have a linear dependence on Fe/Mg component. Both 1 µm and 2 µm absorption bands seem to follow approximation lines in the three-dimensional spaces of center wavelengths, Ca and Fe components. The successful application of the exchange Monte Carlo method to a wide range of clinopyroxenes would have a potential to expand the applicability of MGM to a variety of space/ground-based observations, especially when we cannot rely on prior information of the mineralogy.
where C 0 and C 1 are respectively intercept and slope of the continuum in frequency space.
The technical difficulty of applying the modified Gaussian model (MGM) is due to a local minima problem. The numerical algorithm of the widely used MGM utilizes the steepest descent method [36,38] or total inversion algorithm [37,39], both of which are not guaranteed to converge into a global solution. The analyzer can find an optimal solution only when the appropriate initial parameters are provided, based on preliminary knowledge of mineralogy [40,41]. With inaccurate initial parameters, however, these gradient descent methods converge into local minima, and thus the analyzer must manually adjust initial parameters and calculate the model iteratively to obtain the desired solution [16,42]. This would be a significant obstacle when one needs to automatically analyze large spectral databases obtained by space missions. In addition, preliminary knowledge of mineralogy may not always be available for space/ground-based observations, especially when the reflectance spectra are the only useful obtained data for interpreting the mineralogy of target bodies. Given the recent rapid increase of reflectance spectral data, automation of deconvolution analysis without requiring preliminary information on mineralogy is warranted. Makarewicz et al. [43] and Parente et al. [42] recently developed an algorithm to select initial band parameters automatically, based on inflection points of the derivatives of observed spectra. Although their algorithm does not depend on prior information, many spurious local minima and inflection points due to noise lead the authors to apply a smoothing filter, yielding arbitrariness on their analyses.
In order to overcome the local minima problem, the exchange Monte Carlo method, also known as parallel tempering [44], has been widely applied in the fields of physics, chemistry, biology, engineering and materials science [45]. Nagata et al. [46] developed a Bayesian spectral deconvolution model combined with the exchange Monte Carlo method with application to visible to near-infrared (Vis/NIR) reflectance spectra of fayalite and forsterite. This method is an improved algorithm of the Markov chain Monte Carlo method, aimed to avoid local minima traps [47] and to remove the arbitrariness originated from initial parameters. In order to solve the local minima problem, the simulated annealing scheme [48] has been incorporated into the model of Nagata et al. [46]. This algorithm introduces a pseudo-temperature and attempts to find the global minimum by heating and cooling the system. Nagata et al. [46] showed that the method can deconvolve reflectance spectral data of fayalite and forsterite into a few Gaussians with a continuum, purely based on the observed data, without requiring preliminary information of the band structure of olivine absorptions. In this paper, we report the applicability of the exchange Monte Carlo method to more complex rock-forming minerals (i.e., clinopyroxene). As described below, since the behavior of pyroxene spectra with the change of chemical composition is relatively well understood, the use of reflectance spectra of pyroxene minerals is suitable for testing the new spectral deconvolution method. Clinopyroxene (Cpx), with its general formula being (M2)(M1) (SiAl) 2 O 6 , is one of the most important rock-forming mineral groups due to both its rich abundance on solid bodies in the solar system and distinguished absorption features [20,34,49]. Cpx includes a wide range solid solution of Mg, Fe and Ca compositions and has two crystal structures of C2/c and P2 1 /c [50], which could reflect various physical and chemical processes inside planetary bodies, such as the thermal history of magma [51]. Due to its wide range of chemical compositions, reflectance spectra of Cpx minerals vary significantly, and are generally grouped into three types: type-A, B and A/B [20,52,53]. Three major absorption bands are observed in type-B spectra, centered around 1.0, 1.2, and 2 µm, attributed to spin-allowed crystal field transitions of Fe cations in the octahedral (M1 and M2) sites [20]. These band centers are known to vary due to total iron and calcium content [19,[52][53][54][55]. On the other hand, type-A spectra lack a strong 2 µm band, interpreted as a low Fe 2+ content in the M2 site [53]. Type-A/B spectra are intermediate between type-A and B, although the boundaries are not well defined. MGM analyses have been performed to natural Cpx [36,53], synthetic Cpx [56], and mixtures of orthopyroxene-clinopyroxene [37,41].

Exchange Monte Carlo method
The technical details of the exchange Monte Carlo method is described elsewhere [46], thus we briefly summarize the key parameters of the model. In our study, the hyperparameters for Gamma and Gauss distributions used to yield probability densities of the parameters are those identified in Nagata et al. [46]: η a = 3.0, λ a = 2.0, ν 0 = 1.25, ξ 0 = 2.5, η b = 5.0 and λ b = 0.04. The total number of temperatures L in our study was 80, and the inverse temperature β l given by: Figure 1 shows typical examples of the evolution of root mean square (RMS) during the exchange Monte Carlo calculations for sample 088 using 4 Gaussian functions. Rapid decreases of RMS for the lowest temperature can be observed at about Monte Carlo steps = 300, 500, 700, 1500, and 6000, due to the parameter exchange between the lowest and middle temperatures. From Figure 1, it can be observed that the RMSs of higher temperatures are generally larger than those of lower temperatures, showing attempts to find better global minimum with wider fluctuations. On the other hand, for lower temperatures, each iteration attempts to find local minimum within a parameter range narrower than higher temperatures. After the Monte Carlo step exceeds 10 4 , although the model still attempts to find better solution Gaussian functions. Significant reduction of RMS for the lowest temperature can be observed approximately at Monte Carlo steps: 300, 500, 700, 1500, and 6000, which shows that the parameter exchange between the lowest and middle temperatures have occurred. The first 100,000 steps were used for the burn-in period and the last 20,000 steps for the expectation value calculations.
with certain fluctuations of parameters, the RMSs remain almost constant, and the expectation values of parameters converge to the best solution. Similar to Nagata et al. [46], the Monte Carlo calculation was iterated through 100,000 steps for the burn-in period and 20,000 steps for the expectation value calculations. Errors of band parameters are estimated from 2σ based on ten runs using a different series of random numbers.
The number of Gaussian functions, K, is an important parameter for deconvolution analysis. The model usually improves with more Gaussians. Though, too many Gaussians may cause overfitting, and the solution can be physically unrealistic. We performed spectral deconvolution using a various number of Gaussian functions. For one spectrum, we varied K from 3 to 10, thus the total number of free parameters ranges from 11 to 32, including the intercept and slope of the continuum. In order to select an optimal K for the deconvolution, we calculate the free energy, or stochastic complexity [57,58], which is an evaluation function for the model section problem [46]. With increasing K, we find that the free energy decreases, and when K is higher than about 5, it converges and fluctuates around a low value. Thus, we chose minimum Ks from the region where the free energies converge and interpret them as optimal Ks, listed in Table 1.

Spectral data
We used the visible to near-infrared spectra of synthetic clinopyroxene samples with a wide compositional range collected at the KECK/NASA Reflectance Experiment Laboratory (RELAB) at Brown University [59,60]. The sample IDs for the RELAB catalogue is summarized in Table 1. In order to compare our results with conventional MGM analyses, we collected 31 reflectance spectra that have been analyzed by a previous study [56]. The method of synthesis is detailed in Turnock et al. [61]. Individual synthetic pyroxene grains typically have 15-25 µm in size, however, these grains formed clumps [56]. Thus, samples are crushed and sieved at <45 µm. The spectra were measured at 5 nm intervals over the wavelength range of 0.3-2.6 µm. The incidence and emission angles were 30° and 0°, respectively [56,59,60]. We performed spectral deconvolution only over the wavelength range of 0.4-2.6 µm, since the standard deviations become larger near the shorter wavelengths [46]. The chemical compositions of the samples are measured with electron microprobe by Klima et al. [56]. The compositions of individual pyroxene samples are indicated in molar ratio with endmember compositions of enstatite (En: Mg 2 Si 2 O 6 ), ferrosilite (Fs: Fe 2 Si 2 O 6 ) and wollastonite (Wo: Ca 2 Si 2 O 6 ) and plot on a pyroxene quadrilateral (Table 1 and Figure 2). Minor compositions typical for natural CPx samples, such as Cr, Mn, Al and Fe 3+ , are not observed [56]. Only the mineral structure of sample 088 is reported to be P2 1 /c [62], while the mineral structures for the remaining samples are not available. Based on the nomenclature of clinopyroxene [50], we assumed low-Ca pyroxene specimens with Wo < 20 to be pigeonite, high-Ca pyroxene with Wo > 45 and En > 25 diopside, and high-Ca   and pressure for a geological timescale [63]. All of the reflectance spectra of synthetic pyroxene samples are shown in Figure 3. Spectra of pigeonite and augite are shown in Figures 3A-3C. The band minima of 1 μm absorptions for pigeonite and augite locate slightly shorter than 1 μm, similar to synthetic orthopyroxene [64]. Their band minima of 2 μm absorptions, caused by spin-allowed crystal field transition of Fe 2+ in the M2 site [52,65], locate at longer than 2 μm similar to Ferich low-Ca pyroxene (orthopyroxene). Pigeonite and augite also show distinctive 1.2 μm absorption bands, which is attributed to spinallowed crystal field transition of Fe 2+ in the M1 site [66]. In Figures 3A-3C, the spectra of pigeonite and augite specimens, which are assigned to type-B spectra [52], show distinctive 1 and 2 μm absorption bands. On the other hand, some of the high-Ca pyroxene (i.e., diopside and hedenbergite) lack distinctive 2 μm absorption bands, as shown in Figures 3D and 3E. They are assigned to type-A pyroxene, in which the M2 sites are saturated cations other than Fe 2+ , such as Ca 2+ [52]. The spectrum of sample 083 has a broad 1 μm absorption band, probably due to a composite absorption by M1 bands near 1 and 1.2 μm [56]. This type of spectrum, exemplified by sample 083, may not be a suitable subject for MGM analyses, since the shape of 1 μm absorption band is far from Gaussian. However, since our objective in this paper is to test the applicability of the exchange Monte Carlo method and to compare the results with those of conventional MGM analyses, we included sample 083 in our analysis. Fitting additional small absorption bands centered near 0.50 and 0.55 μm in type-B spectra, which are associated with spin-forbidden crystal field transitions [66], is beyond the scope of this paper, since our focus was to evaluate the applicability of the model to two major absorption bands of clinopyroxene centered about 1 and 2 μm.

Deconvolution with the exchange Monte Carlo method
Deconvolution results of pyroxene spectra are shown in Figure 4. The best optimized parameters for 1, 1.2, and 2 μm bands are summarized in Table 1. Errors estimated from 2σ based on ten runs using different series of random numbers are also shown. Most of the spectra are best fitted by 4 to 6 Gaussians with appropriate continuum. For example, sample 053 is fitted with 4 Gaussians with their centers locating at 0.091, 0.965, 1.258, and 2.180 μm and a nearly constant continuum. For pigeonite and augite, our results show that the spectra are mainly fitted with 3 Gaussians with their centers locating at about 1, 1.2, and 2 μm, and with one Gaussian in the wavelength region ranging from ultraviolet to visible and with a continuum being nearly constant or gradually decreasing toward shorter wavelengths (Figures 5a-5d). Each Gaussian function could be physically interpreted to represent spin-allowed crystal field transitions with 1 μm for Fe 2+ in the M1 and M2 sites, 1.2 μm for Fe 2+ in the M1 site and 2 μm for Fe 2+ in the M2 site (Figures 6a-6d). A Gaussian in the wavelength region ranging from ultraviolet to visible could represent oxygen-metal charge transfers which are centered within the ultraviolet region [34] (Figure  7). Although additional small bands centered near 0.7 μm may be necessary to improve the fitting such as shown in sample 054, we find that the additional band does not significantly affect the band center for 1, 1.2, and 2 μm absorptions. Our deconvolution results for pigeonite and augite are consistent with those by Klima et al. [56].
Deconvolution analyses for type-A spectra including some of the diopside and hedenbergite are not as straightforward as type B, since the 1 μm bands are too narrow or too wide for a single Gaussian to fit, and some of the 1 and 2 μm bands are asymmetrical. The non-Gaussian shape of type-A spectra results in larger errors of band parameters when compared to type-B spectra (Table 1). Also, more Gaussian functions are needed to fit type-A spectra, reflecting the complex shape of the spectra. For example, samples 037 and 082 were fitted with 7 Gaussians, while sample 083 with 10 Gaussians of which 4 Gaussians are centered below 0.5 μm. Comparing our fitting results of sample 082 and 083 with those reported by Klima et al. [56], we find that the exchange Monte Carlo method yields more symmetrical configurations for Gaussians in 1 μm band. We also find that our optimal deconvolution results for samples 082 and 083 do not require the very weak absorption bands which were assigned as M2 absorptions in Klima et al. [56]. These M2 absorptions in Klima et al. [56] are significantly weak compared with the strongest M1 absorption in 1 μm band, and the M2 absorption near 1 μm was mostly covered with the M2 absorptions. For type-A spectra, it is difficult to resolve the discrepancy between the results obtained by Klima et al. [56] and this study, however, since both of the modeling results can reproduce the observed spectra almost equally well. Nevertheless, the large errors of band parameters for type-A spectra indicate that caution should be taken when performing spectral deconvolution for type-A spectra using MGM. Thus, being able to estimate errors in the fitting results based on random initial parameters is also an advantage of the exchange Monte Carlo method for assessing the statistical robustness of fitting results, which has not been performed by previous MGM analyses.

Band shift as functions of Ca, Fe, and Mg Contents
The center positions of Gaussian functions corresponding to 1, 1.2, and 2 μm absorptions are summarized in Table 1. Position shifts of 1  Figure 5, while position shifts of 2 μm band are shown in Figure 6. Our results are found to be consistent with those of Klima et al. [56], suggesting that the exchange Monte Carlo method developed in this study is able to extract results similar to those obtained by conventional MGM analysis.
The band centers of 1 μm absorption have a linear dependence on Wo (Ca/(Mg+Fe+Ca) molar%) contents, moving to longer wavelengths with increasing Wo. For diopside and hedenbergite with Wo ≈ 50, the band centers of 1 μm absorption largely scatter within the range of 1.03-1.08 μm. Pigeonite also shows a large variance as a function of Wo component, however, the band centers of 1 μm absorptions have a linear dependence on En (Mg/(Mg+Fe+Ca) molar%) or Fs (Fe/(Mg+Fe+Ca) molar%) components. The band centers of 1 μm absorption of augite scatter widely as a function of Fs content, and distinct dependence on Fs content were not observed. For diopside and hedenbergite samples, with the exception of sample 043, the band centers of 1 μm absorptions seem to have a linear dependence on Fs content, moving to longer wavelengths with increasing Fs content. En content seems to have no obvious influence on the band center of 1 μm absorption, except for pigeonite. We find that pigeonite, augite and diopside-hedenbergite can be separated with the use of 1 μm band position as a function of Fs or En content, as shown in Figure 5.
The band centers of 2 μm absorption seem to follow an approximation line on the space of Fs and Wo contents (Figure 6a). The band centers of 2 μm absorption for pigeonite and augite move to longer wavelengths with an increase of Wo content. The average center wavelengths of 2 μm absorptions for diopside and hedenbergite remain almost constant at around 2.3 μm, but they scatter significantly around the average value, reflecting the asymmetry shape of 2 μm absorptions for some of the type-A spectra. The band centers of 2 μm absorptions of pigeonite move longer wavelengths with increasing Fs content. Overall, the band center of 2 μm absorption is separated between low-Ca clinopyroxene (pigeonite) and high-Ca pyroxene (augite, diopside, and hedenbergite) with a gap from 2.20-2.25 μm. Klima et al. [56] interpreted the gap to be a transition zone of mineral structure between P2 1 /c and C2/c. Our results generally agree with their interpretation with one exception, i.e., sample 056 is assumed to be pigeonite in Klima et al. [56].
The diagram between 1 µm and 2 µm band positions is shown in Figure 7. Generally low-Ca pyroxene locates in a shorter wavelength region, while high-Ca pyroxene locates in a longer wavelength region [19,41,52,56]. Our results are consistent with previous studies based on MGM analysis [41,56].

Discussion
In order to avoid the local minimum problem, a Bayesian spectral deconvolution method with the exchange Monte Carlo algorithm has been applied to visible to near infrared reflectance spectra of synthetic Cpx with wide ranging Mg, Fe, and Ca contents. The results obtained in this study generally agree well with conventional MGM analyses. Here, we discuss some potential interpretations for the deconvolution results, as well as the discrepancies with previous results obtained by MGM analyses.

Fitting type-A spectra
Type-A spectra generally appear in high-Ca pyroxene such as diopside and hedenbergite, although previous studies suggest that there is no simple relationship between chemical composition of natural pyroxenes and type-A spectra [20,52,53]. The spectral data of samples 076, 037, 083, and 036 are categorized in type-A, which lack distinctive 2 µm bands. Since the absorption near 2 µm is mainly caused by spinallowed crystal field transition of Fe 2+ in the M2 site, the absence of 2 µm band from type-A spectra is interpreted as a result of the replacement of Fe 2+ in M2 site by other cations such as Ca 2+ and Fe 3+ [52,67]. On the other hand, spin-allowed crystal field transition due to the M1 site, with the absorption center locating near 1 and 1.2 µm, appear strong, and the two bands are combined to yield single broad band around 1 µm. It is difficult to fit type-A spectra with manually provided initial parameters (e.g., center, strength of Gaussian) using a conventional MGM algorithm, since the final solutions directly depend on the initial parameters, especially for the broad 1 µm band due to its non-Gaussian shape [56]. Although it is still difficult for the exchange Monte Carlo method to delineate the optimal number of Gaussians, we find that the center wavelengths of Gaussian functions for the broad 1 µm band do not move significantly with varying K. For example, the broad 1 µm band of sample 083 was fitted with two major Gaussian functions which could correspond to crystal field absorptions due to the M1 site whose center wavelengths of absorptions locate near 1 and 1.2 µm. Although the estimated errors are large, the overall result is consistent with previous deconvolution analyses of type-A spectra [53,56], suggesting that even with manually provided initial parameters, previous analyses obtained statistically optimal results. Since manually fitting type-A spectra is more difficult than type-B spectra, the exchange Monte Carlo method can be a useful tool for future deconvolution analysis for type-A spectra, of which are likely to appear in high-Ca pyroxene [52,53]. In addition, natural Cpx incorporate many minor compositions such as Al, Ti, Mn and Cr [49,50]. Such minor elements in natural Cpx yield more complex spectra than synthetic Cpx, thus the scheme presented in this paper can be useful to deconvolve such spectra. Figure 5a indicates that the Cpx seem to follow a linear relationship among the 1 µm band position, with Fs and Wo contents. As shown as broken lines in Figures 5b and 5d, linear relationships are observed for pigeonite and diopside-hedenbergite with their Fs or En content. It should be noted that for diopside-hedenbergite, linear approximation was performed for all of the samples except sample 043. Sample 043 was omitted as an outlier because the spectrum shows higher reflectance and weaker absorptions compared with other samples, generally suggestive of the effect of glass [16,56,68]. Although it has been well documented for low-Ca pyroxene that the linear dependence of 1 µm band centers on Fs or En content [19,41,52,54,56], the relationship between the 1 µm band center and Fs or En content for high-Ca pyroxene has been poorly constrained by MGM analyses. Clénet et al. [41] analyzed only two high-Ca natural pyroxene samples in the diopside-hedenbergite region, thus no clear relationship has been derived. Although Klima et al. [56] analyzed 13 high-Ca pyroxene samples in the diopsidehedenbergite region; they observed no clear relationship between the 1 µm band center and Fs or En content. Our analyses show that the 1 µm band center of high-Ca pyroxene (i.e., diopside and hedenbergite) depends on Fs or En content with the band center moving longer wavelengths with increasing Fs content. This tendency is similar to the relationship between the 1 µm band center and the Fe/Mg ratio of Cafree orthopyroxene [64]. For Ca-saturated synthetic Cpx in which most of the M2 site is dominated by Ca 2+ , varying Fe/Mg ratio would affect only Fe/Mg in the M1 site because the M2 site is dominated with larger Ca 2+ cations. Thus, the change of Fe/Mg would appear only in the 1 µm band center caused by a crystal field transition in the M1 site, but not in the 2 µm band center caused by the absorption due to the M2 site. We note, however, that large errors are included in the modeling results for type-A spectra (Table 1), thus additional data is necessary to confirm our interpretation.

Band center in three-dimensional spaces
As mentioned above, both the 1 µm and 2 µm bands seem to plot on approximate lines in the three-dimensional (3D) space of Fs-Wo content (Figures 5a and 6a). Low-Ca pyroxene samples (pigeonite) used in this study locates in the high-Fe and low-Ca regions with both the band centers being at shorter wavelengths. On the other hand, high-Ca pyroxene samples (diopside and hedenbergite) used in this study locates in the low-Fe and high-Ca regions with both the band centers being at longer wavelengths. Augite with intermediate-Wo contents distribute a in a linear trend between low-Ca and high-Ca pyroxene. Each mineral group displays three distinct clusters on these 3D spaces. The linear trend seen in both 1 µm and 2 µm bands could be understood from a crystallographic point of view. Whereas high-Fe and low-Ca Cpx samples contain more Fe 2+ cations in both the M1 and M2 sites, leading to a decrease of the bond lengths due to the smaller size of Fe 2+ , which results in an increase in the crystal field splitting [34,69,70], low-Fe and high-Ca Cpx samples have larger Ca 2+ cations, which dominate the M1 and M2 sites, leading to an increase in the bond lengths and a decrease in the crystal field splitting.

Application to future remote reflectance spectroscopy
Despite requiring meticulous parameter adjustment and prior information of mineralogy, MGMs using gradient descent method have been applied to reflectance spectra not only of laboratory data but also remote sensing data of the Moon [71], Mars [4,6,[72][73][74][75], and some asteroids [76,77]. By applying the exchange Monte Carlo method to spectra of synthetic Cpx samples, we have shown that it is able to yield results consistent with both conventional gradient descent methods and crystal field theory. This means that the exchange Monte Carlo method could be applicable to at least some of the previous remote sensing data which have been analyzed using conventional MGM methods. Because the application of MGM analyses has been limited due to meticulous parameter adjustment, the exchange Monte Carlo method may have a potential to expand the applicability of MGM to a variety of space/ground-based observations, especially when we cannot obtain prior information of the target body beforehand. Remote sensing data of terrestrial reflectance spectroscopy have been validated with various reference spectra observed by in situ sample analyses, increasing the confidence of interpretation of remote sensing data [12,17,78,79]. However, such an analysis is not always available due to various geological contexts. The situation is more serious for space/ ground-based observations of small bodies in the solar system, because reflectance spectra are the only available compositional information for most of the small bodies such as Phobos, Deimos [80,81] and asteroid 1999 JU3 which is the target of Japanese space mission Hayabusa 2 [82]. In fact, asteroids are clustered based solely on reflectance spectra, although their detailed mineral compositions are poorly constrained [83]. Considering the large number of small bodies in the solar system, it is unlikely to send probes to each small body to conduct in situ sample analysis or sample return. Given the limited data and resource, interpreting reflectance spectra without assuming a priori information is essential not only for planetary science but also for mission planning. The exchange Monte Carlo method could be significantly useful under such circumstances, where meaningful compositional information other than reflectance spectra is not available.

Conclusions
We applied a Bayesian spectral deconvolution method with the exchange Monte Carlo algorithm to visible/near infrared reflectance spectra of synthetic clinopyroxene of diverse compositional variation, in order to avoid the local minimum problem and to remove the arbitrariness originated from initial parameters, inherent in the previous Modified Gaussian model. Our results indicate that the exchange Montel Carlo method is able to yield the consistent results obtained by conventional Modified Gaussian model. Since our model does not rely on a preliminary knowledge of the reflectance spectrum of mineral, the results suggest that previous spectrum analyses have obtained statistically optimal results. The successful application of our model to reflectance spectra of minerals indicates that this model could be applied to an automatic deconvolution analysis for a large spectral database, especially useful for space missions when preliminary knowledge of mineralogy is not available. Given the recent and future advancement of space missions, the exchange Monte Carlo method could be a useful tool for analyzing a wide range of minerals and remote sensing data of rocky bodies in the solar system.