{'fraction_non_alphanumeric': 0.03512232673386418, 'fraction_numerical': 0.024056388173879575, 'mean_word_length': 4.466140696909927, 'pattern_counts': {'":': 0, '<': 4, '<?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': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Non-parametric Approach to Constrain the Transfer Function in Reverberation Mapping

Broad emission lines of active galactic nuclei stem from a spatially extended region (broad-line region, BLR) that is composed of discrete clouds and photoionized by the central ionizing continuum. The temporal behaviors of these emission lines are blurred echoes of the continuum variations (i.e., reverberation mapping, RM) and directly reflect the structures and kinematic information of BLRs through the so-called transfer function (also known as the velocity-delay map). Based on the previous works of Rybicki&Press (1992) and Zu et al. (2011), we develop an extended, non-parametric approach to determine the transfer function for RM data, in which the transfer function is expressed as a sum of a family of relatively displaced Gaussian response functions. Therefore, arbitrary shapes of transfer functions associated with complicated BLR geometry can be seamlessly included, enabling us to relax the presumption of a specified transfer function frequently adopted in previous studies and to let it be determined by observation data. We formulate our approach in a previously well-established framework that incorporates the statistical modeling of the continuum variations as a damped random walk process and takes into account the long-term secular variations which are irrelevant to RM signals. The Application to RM data shows the fidelity of our approach.


INTRODUCTION
The broad-line regions (BLRs) in active galactic nuclei (AGNs), with sizes of light days to light weeks (Kaspi et al. 2000;Bentz et al. 2013), are too compact to be spatially resolved with the existing telescopes, except in a few cases in which the BLRs are magnified by gravitational lensing (Sluse et al. 2012;Guerras et al. 2013). The reverberation mapping (RM) technique provides a powerful and unique way to probe the structures and kinematics of BLRs through analyzing the temporal variation patterns of emission lines and the ionizing continuum (Blandford & McKee 1982;Peterson 1993). Over the past four decades, RM studies have successfully measured BLR sizes for ∼ 60 nearby Seyfert galaxies and quasars (e.g., Bentz et al. 2013;Du et al. 2015). Particularly influential is the discovery that BLR sizes closely correlate with optical luminosities of (sub-Eddington) AGNs, as expected from simple photoionization theory (e.g., Kaspi et al. 2000Kaspi et al. , 2005Bentz et al. 2013;Du et al. 2015). This relationship, in combination with the velocity widths of broad emission lines, opens up an economic way for the virial mass estimation of supermassive black holes (e.g., Vestergaard & Peterson 2006;Ho & Kim 2015). Meanwhile, velocity-binned RM analysis for several AGNs with highquality spectroscopic data has crudely seen signatures for diverse kinematics in BLRs, including radial inflows/outflows and non-radial, virialized Keplerian rotations Denney et al. 2010;Grier et al. 2013b;Peterson 2014;Du et al. 2016).
Despite these advances, the detailed information for structures and kinematics of BLRs remain largely unknown. The major challenges are from RM data sets reliably inferring the transfer function or velocity-delay map and from obtaining the structure and kinematic information of BLRs from this map. Traditional cross-correlation analysis (e.g., Bentz et al. 2013) or the previous method using a simple specified transfer function (e.g., Rybicki & Kleyna 1994;Zu et al. 2011Zu et al. , 2016 yields no more than the characterized BLR sizes. Great effort has gone into developing more sophisticated methods. There are generally two classes of methods depending on the strategy adopted in solving the integral equation that links the variations of continuum and emission lines in RM. One class relies on directly inverting the integral equation, including early attempts based on the Fourier transform (Blandford & McKee 1982) and regularized linear inversion (RLI) methods (Krolik & Done 1995;Skielboe et al. 2015). The method using the Fourier transform was found to require high-fidelity data and thus has not seen much application (Maoz et al. 1991), whereas RLI seeks to discretize the integral equation and employ a differencing operator to suppress the noise. The downside of this class of methods is that the measurement noises cannot be self-consistently incorporated, thus, this class may be very noise sensitive in some cases. The other class relies on indirectly inferring the "best" transfer function that fits the data by supposing prior limits or models on it. The maximum entropy method (Horne 1994) and BLR dynamical modeling method (Pancoast et al. , 2014aLi et al. 2013) belong to this class. The maximum entropy method finds the "simplest" solutions for the transfer function by maximizing the entropy on the premise of a reasonable fit to the data. However, it does not allow for straightforward uncertainty estimates or model selection (Horne et al. 2003). The BLR dynamical modeling method presumes a specific but sufficiently flexible toy model for BLRs and directly infers the model parameters; therefore, the transfer function is only an intermediate product FIG. 1.-Schematic of the transfer function for a system that consists of discrete clouds. The transfer function can be expressed as a sum of many relatively displaced response functions of local adjacent clouds. For illustrative purposes, the red (blue) curve in the transfer function corresponds to the response from the clouds highlighted in red (blue).
-Illustrative examples for three types of BLR geometry in which the transfer functions are approximated by a family of Gaussian response functions: (left) a spherical shell, (middle) two disk rings, and (right) an inclined disk ring. Arrows indicate the line of sight. In the bottom panels, black lines represent the realistic transfer functions, and blue lines represent the sum of the Gaussian functions (gray dotted lines). Our approach can particularly cope with the cases of multimodal or asymmetric transfer functions (e.g., middle panels), in which the previous JAVELIN method may encounter difficulties (see Appendix C). Li et al. 2013). It cannot yet satisfactorily assess the systematic errors of the adopted model.
In this paper, we describe an extended, non-parametric method to determine the transfer function from RM data, based on the previous works of Rybicki & Press (1992) and Zu et al. (2011). Motivated by the scenario that BLRs are composed of a large number of clouds so that a local adjacent collection of these clouds effectively respond to the central ionizing continuum as a Gaussian response function, we express the transfer function as a sum of a family of relatively displaced Gaussian functions. Therefore, we can account for arbitrary shapes of the transfer function, thereby largely relaxing the need for presuming the transfer function (e.g., Rybicki & Kleyna 1994;Zu et al. 2011) or BLR models (e.g., Pancoast et al. 2011;Li et al. 2013). In this respect, our method is model independent and is particularly able to cope with complicated shapes of transfer functions (e.g., multimodal) for diverse BLR structures and kinematics.
The paper is organized as follows. We describe the methodology of our approach in detail in Section 2 and perform simulation tests in Section 3. Section 4 presents application of our approach to three RM AGNs to illustrate its fidelity. A discussion and conclusions are given in Section 5. All the time delays quoted in the text are given in the observer's frame.

The Transfer Function
In RM, variations of broad emission lines f l (t) are blurred echoes of continuum variations f c (t) through the transfer function Ψ(τ ) as (Blandford & McKee 1982;Peterson 1993) where the integral limits are set to account for all possible positively and negatively delayed responses (e.g, Marshall et al. 2008;Shappee et al. 2014). Here, we show only velocityunresolved RM so that velocity information is not included in the transfer function. However, it is easy to extend our approach to velocity-resolved RM by applying Equation (1) to different velocity bins of the emission lines (see below for a discussion). In addition, we also neglect the nonlinearity of the response for simplicity (Li et al. 2013). Nevertheless, we will provide below a way to incorporate non-linear response. Mathematically speaking, it is known that an arbitrary function can be expressed as a sum of many slightly displaced is the Dirac delta function. In practical calculations, one instead uses a set of basic functions with spatial extent to describe a smooth function. A general adoption is Gaussian functions, which are often found to be adequate (Sivia & Skilling 2006, p.140). In this spirit, the transfer function in RM is now written as where ω represents the common width of Gaussian responses, τ k and f k are the mean lag and weight of kth response, and K is the number of responses, which can be regarded as a "smoothing parameter." The Gaussian functions are fixed to a common width for the sake of simplicity. Such a choice also helps to relax the parameter degeneracy and improve the computational efficiency of the Bayesian inference described below. The physical interpretation for such an expression is straightforward. The system (i.e., BLRs) under consideration is composed of a large number of clouds. Each cloud responds to the variation of the central source as a δ−response function. A collection of these local adjacent clouds effectively produce a Gaussian response function with a width. This is similar to the principle of the local optimally emitting cloud model, widely applied in photon-ionization calculations (Baldwin et al. 1995). Figure 1 shows a schematic illustration of the transfer function of BLRs. In Figure 2, we show examples for three types of BLR geometry-a spherical shell, double disk rings, and an inclined disk ring-in which the transfer function is approximated by a family of Gaussian functions. In particular, there is evidence that Hβ BLRs have two distinct components (e.g., Brotherton 1996;Sulentic et al. 2000;Hu et al. 2008, 2012 andreferences therein), which indicate that their transfer functions may be bimodal, as illustrated in the middle panels of Figure 2.

Bayesian Inference
With Equation (2), we can directly use the framework developed by Rybicki & Press (1992) and Zu et al. (2011) to infer the best transfer function. We stress that all the formulae presented in this section have been well established by Zu et al. (2011). We list them below for the sake of completeness.
Continuum variation can be well described by a damped random walk (DRW) process (e.g., Kelly et al. 2009Kelly et al. , 2014Zu et al. 2011;Li et al. 2013). A DRW process is a stationary process such that its covariance function at any two times t i and t j only depend on the time difference ∆t = t i − t j , which can be simply prescribed as where angle brackets represent the statistical ensemble average, f c (t) is the continuum driven by the DRW process, τ d represents the typical damping timescale, and σ d represents the standard deviation of the process on a long timescale (≫ τ d ).
On a short timescale, the variation amplitude of the process is σ d 2t/τ d . If this is combined with Equation (1), the covariance between the line and the continuum can be written as and the covariance of the line is In Appendices A and B, we present a detailed derivation for the analytical expressions of these covariances in terms of the error function. The Bayesian posterior probability is obtained as follows (see Press 1992 andZu et al. 2011 for a thorough derivation). Let y be a column vector comprised of the light curves of both the continuum and the line. A set of measurements in a campaign is equal to the realization of some underlying signals with measurement noises as where s is the signal of the variations, of which the portion for the continuum is described by a DRW process, and n represents the measurement noises. Here the term Lq is usually used to model general trends linearly varying with time in the light curves (Rybicki & Press 1992), where L is a matrix of known coefficients and q is a column vector of unknown linear fitting parameters. For example, to remove separate means of the light curves with a total of N + M measurements, one just configures L to be a 2 × (N + M) matrix, where N and M are the number of measurements for the continuum and emission line, respectively. The matrix L has entries of (1, 0) for the continuum data points and (0, 1) for the line data points (Rybicki & Press 1992;Zu et al. 2013). If one desires to detrend the light curves with long-term secular variations using one-order polynomials, one needs to configure L as a 4 × (N + M) matrix (see below for details). Unless stated otherwise, we remove the means of the light curves by default in the following calculations. The signal s and noise n are assumed to be Gaussian and independent, so that their probabilities are simply written as and respectively, where S = ss is the covariance matrix of s given by Equations (3)-(5) and N = nn is the covariance matrix of n. There is no prior information for the linear fitting parameters q, and we assign the vector a uniform prior probability. Now the probability for a realization y is where θ represents the parameter set for the transfer function in Equation (2) and the probability P(q) is uniform, so that it is neglected in the above integral. With some mathematical implementation, it can be shown that (Zu et al. 2011) where C = S + N and C −1 and the best estimate for the linear fitting parameter iŝ The corresponding variance of the best estimate is, respectively, and The best estimate for the light curves iŝ and the corresponding variance is Note that Equations (11), (12), and (15) also provide a way to estimate the value of the light curves at some specified time that may not be one of the already measured points. According to the Bayes's theorem, the posterior probability distribution for the free parameter set is given by where P(y) is the marginal likelihood that serves as a normalization factor to the posterior probability and P(σ, τ d , θ) is the prior probability for the parameters. We maximize this posterior distribution to obtain the best estimate for the free parameters using a Markov-chain Monte Carlo (MCMC) method as described in the next section. Long-term secular variability is incidentally detected over the duration of the campaign in some RM objects (e.g., Denney et al. 2010;Li et al. 2013;Peterson et al. 2014). Such long-term trends are uncorrelated with reverberation variations and may bias the reverberation analysis (Welsh 1999). The linear fitting parameter q provides a natural way to remove such trends. If the long-term trends vary linearly with time (q 1 + q 2 t), the matrix L has the form of (Rybicki & Press 1992;Zu et al. 2011) where t c i is the ith time of the N measurement points of the continuum and t l i is the ith time of the M measurement points of the line. In this case, q is a column vector with four entries, which represent the coefficients of the linear trends for the continuum and line data. For the long-term trends with high-order polynomial variation, one can similarly write out the matrix L by just adding extra columns. In our calculations, when the light curves of the emission line and the driven continuum visually show different long-term trends, we switch on the detrending procedure by assuming the trends vary linearly with time (see below for application to the RM data of Mrk 817).

Markov-Chain Monte Carlo Implementation
The parameters in our approach include two parameters (σ and τ d ) pertaining to the DRW process in Equation (3) and 2K + 1 parameters ( f k s, τ k s, and ω) pertainng to the transfer function in Equation (2). To relax parameter degeneracy, we introduce the following simplifications.
• The central lags τ k of Gaussian responses are assigned by a uniformly spaced grid of time lag over a range of interests, whereas the amplitude f k is free to determine. Specifically, we set τ k = τ 1 + (k − 1)∆τ for 1 k K, where ∆τ is the grid space and (τ 1 , τ K ) is the range of the time lag. Given with the range of the time lag, the value of ∆τ is determined by the number (K) of Gaussian functions, i.e., ∆τ = (τ K − τ 1 )/(K − 1). The degree of freedom is now reduced to K + 3.
• The width of Gaussian responses ω is limited by two facts. If ω is too small, unnecessary structures in transfer function will become admissible, and hence, the reconstructed light curves will appear to be unrealistically noisy; if ω is too large (e.g., ω > ∆τ ), the neighboring Gaussian responses are indeed indistinguishable. In practice, we find it is adequate to adopt a prior limit of ω ∼ (∆τ /2, ∆τ ).
The prior probability P(σ, τ d , θ) in Equation (17) is assigned by assuming that all the free parameters are independent of one another. The Gaussian width has a uniform prior, and the rest of the parameters have a logarithmic prior, which is a usual choice for "scale" parameters that involve a wide range of values (e.g., Sivia & Skilling 2006, p. 108). Lastly, we need to choose the number of Gaussian response functions. Strictly speaking, one should employ the Bayesian model selection to choose the best parameter K, i.e., the one that maximizes the probability P(K|y). However, it is very computationally expensive to calculate P(K|y) since one needs to marginalize all the other free parameters. Therefore, we resort to the relatively more economic "non-Bayesian" approaches based on the maximum-likelihood estimate of the parameters (e.g., Hooten & Hobbs 2015). We use the Akaike information criteria (AIC; Akaike 1973), measures of the relative quality of statistical models. The original form of the AIC is valid only asymptotically; Hurvich & Tsai (1989) proposed a correction to AIC for finite sample sizes (abbreviated to AICc). The AIC is defined by where m is the number of parameters in the model and the hat operator denotes the best estimate for the corresponding parameter. The AICc is defined by where n = N + M is the total number of data points. The best choice of K is the one that minimizes the AICc. We employ an MCMC method to determine the best estimates and uncertainties of free parameters. We use the Metropolis-Hastings algorithm to construct samples of free parameters from the posterior probability distribution. The Markov chain is run in 150,000 steps in total. The best estimates of free parameters are assigned the mean values of  Figure 3. The AICc reache a minimum at K = 13, 10, and 7, respectively. the samples, and the associated uncertainties are assigned the 68.3% confidence levels ("1σ") that enclose the mean values. The uncertainties of the final obtained transfer function in Equation (2) are determined by the error propagation formulae, under the best inferred K. Therefore, the uncertainty in the best inferred K is not included.

SIMULATION TESTS
In order to test our approach, we generate mock light curves of the continuum and broad emission lines with an assumption of the transfer function, and then implement the above procedures to obtain the most probable estimate for the transfer function and compare it with the input one. As outlined by Rybicki & Press (1992), an unconstrained realization of mock light curves is simply a random series that has the covariance matrix C = S + N with specified parameters σ and τ d . We generate random numbers with the covariance matrix C using the Cholesky decomposition (see Zu et al. 2011). The Cholesky decomposition of a matrix is C = MM T , where M is a lower triangular matrix with real and positive diagonal entries. If r is a series of Gaussian random numbers with a zero mean and unity dispersion, then u = Mr has the covariance matrix of C, since uu T = M rr T M T = MM T = C. We first generate a mock continuum light curve and then convolve it with the adopted transfer function to produce a light curve of emission line.

Validity of the Method
We run three simulation tests, presuming the transfer functions to be composed of 1. two displaced top-hats, 2. two moderately mixed Gaussians (so that the transfer function is double-peaked), and 3. two severely mixed Gaussians, respectively. The left panels of Figure 3 plot these three input transfer functions. Specifically, for the first test, the two displaced top hats are (arbitrarily) configured to Ψ(τ ) = 1 15 × 1 for 5 < τ < 10, 2 for 15 < τ < 20, 0 else, where 1/15 is the normalization factor. For the other two tests, the two Gaussians are (arbitrarily) configured to center at 10 and 18 days for the moderately mixed case and 10 and 15 days for the severely mixed case, with the same standard deviation of 2.5 days. The mean time lags for these three transfer functions are 14.2, 14.0, and 12.5 days, respectively, roughly corresponding to a 5100 Å luminosity of ∼ 10 43 erg s −1 according to the relation between the BLR sizes and 5100 Å luminosities of AGNs (Bentz et al. 2013). For the DRW model of all the three tests, the parameter of the damping timescale is set to τ d = 80 days, a typical value for a luminosity of 10 43 erg s −1 (Kelly et al. 2009;Zu et al. 2011;Li et al. 2013); the parameter of the variation amplitude is set to σ d = 1.5 (in arbitrary units) such that the resulting variation amplitude is typically ∼30-50%. The time span of continuum monitoring is 150 days, and the monitoring of the emission line starts 20 days latter. The cadence is typically one day apart. The signal-to-noise ratio (S/N) is set to ∼120. It is worth stressing that the choice of these values is only for illustrative purposes and has no special meaning. The left panels of Figure 3 show the generated mock light curves of the continuum and emission lines for the three simulation tests. Throughout the calculations, the range of the time lag is set to (0-30) days. Since the best number of Gaussians is determined by the AICc, the final transfer function is insensitive to the adopted range of time lag, provided it is broad enough to enclose the time lag of the emission line. The obtained AICc with the number of Gaussian functions K for the three tests are illustrated in Figure 4. The resulting optimal K are 13, 10, and 7, respectively. This means that the optimal K depends not on the sampling cadence but on the shape of the transfer function. The right panels of Figure 3 show the best recovered transfer functions with the optimal K. As can be seen, the structures in the transfer functions are successfully recovered, and all the features in the light curves are also well reconstructed. However, we note that, in the first test, we fail to reproduce the top-hat shape because of the finite sampling and measurement errors of the light curves. Meanwhile, in the second test, although the recovered two peaks are marginally equal, the left peak seems slightly stronger than the right one. These tests suggest that, like the performance of all other RM methods, that of the present approach depends on the data quality and the presence of strong variability features in the light curves; therefore, our approach may be unable to resolve very fine structures of the transfer function in some circumstances.
On the right panels of Figure 3, we also superimpose the interpolation cross-correlation function (CCF) of the light curves using the procedure of Gaskell & Peterson (1987). All the CCFs in the three tests are single-peaked, indicating that CCF analysis hardly reveals more than the characteristic time lags of BLRs. In particular, the CCF in the first test is wide and peaks at τ = 17.5 day. This is an appropriate choice for the time lag if one only concentrates only on a single characteristic lag for the BLR, but clearly a single lag is not sufficient to describe the BLR in this case. In Appendix C, we show the time lag distributions from the JAVELIN method for the above three simulation tests. JAVELIN presumes a top-hat transfer function, and the time lag is assigned the mean lag of the top hat. We illustrate that while the present approach works fairly well in all the three tests, JAVELIN cannot appropriately recover the characteristic lag for the first test with a multimodal transfer function.

Dependence on Data Quality
The above simulations show that our approach can generally well reproduce the input transfer function with optimal data quality. In this section, we test the dependence of the obtained transfer function on data quality with different S/N, sampling cadences, and time spans of the mock light curves. We use for illustration only the double-peaked transfer function-namely, the one used for the second simulation test in the preceding section. We fix the seed for the random number generator so that all the mock light curves in this section have the same variation pattern. This allows us to focus only on the desired factor that influences the obtained transfer function.
In Figure 5, we generate three sets of mock light curves with S/N=120, 80, and 40, respectively. As expected, when S/N is low enough, the weak variability signals in the light curves are overwhelmed by noises, and the information of small structures in the transfer function is missing. Consequently, for the case of S/N=40 (bottom panels of Figure 5), we obtain only a single-peaked transfer function, which has the same mean time lag as the input one.
To check how the sampling rate affects the results, we respectively randomly remove about one-third and two-thirds of the points of the light curves with S/N=120 in Figure 5 and generate two new mock data sets with degraded sampling rates. We show the light curves and the corresponding best recovered transfer function in Figure 6. Again, when the sampling cadence is poorer, more information on the fine structure of the transfer function is lost.
The default time span of the light curves is 150 days. In Figure 7, we reduce the time span to 120 and 90 days, by subtracting the last 30 and 60 days section, of the light curves with S/N=120 in Figure 5, respectively. The mean time lag of the input transfer function is 14 days, so that the overall time spans are still more than six times the mean lag (Welsh 1999). There is no surprise that the transfer function is well reproduced in all the three cases in Figure 7.
In summary, provided the time span is more than several times the mean time lag of the emission line, the S/N and sampling cadence are more critical to recovering the transfer function. We stress again that the overall performance of the present approach depends on the presence of strong variation features in the light curves and the specific structures of the transfer function one desires to resolve. We therefore cannot specify more details about the data requirements for the present approach to faithfully recover the transfer function.

APPLICATION
We select three illustrative RM sources to show the feasibility of our approach: Arp 151, SBS 1116+583A, and Mrk 817. The first two objects were monitored by the Lick AGN Monitoring Project (LAMP 2 ; Bentz et al. 2009;Walsh et al. 2009), and the last object was monitored by Peterson et al. (1998). The RM data of Arp 151 were previously studied in great detail using the maximum entropy method (Bentz et al. 2010), BLR dynamical modeling method Pancoast et al. 2014b), and RLI method (Skielboe et al. 2015). Arp 151 therefore serves as a good example to further verify the validity of our approach. The transfer function of SBS 1116+583A derived by the RLI method of Skielboe et al. interestingly showed a bimodal distribution, meriting an analysis using our approach. The object Mrk 817 shows long-term secular trends in its light curves (Li et al. 2013) and is therefore adopted to illustrate the capability of our approach for coping with detrending.
To compare our results with those of CCF analysis, we define the mean time delay from the transfer function as (Rybicki & Kleyna 1994) The uncertainty of τ m is calculated by the error propagation formulae.

Arp 151
For the sake of comparison with the previous results, we use the B−band photometric light curve of Arp 151 as a surrogate of the continuum as in Pancoast et al. (2014b) and Skielboe et al. (2015). The magnitudes are first converted into fluxes. Since the absolute unit of light curves is no longer important for RM analysis, we adopt an arbitrary zeromagnitude flux for the conversion. We set the range of the time lag for solving the transfer function to be (0-10) days. The best number of Gaussian functions is K = 4 according the AICc.
The top panels of Figure 8 present the reconstruction for the light curves of the continuum and Hβ line and the best recovered transfer function for Arp 151. The mean time lag from Equation (22) is τ m = 3.33 ± 0.26 days, consistent with the previous results from CCF analysis, τ CCF = 4.08 ± 0.60 days ), and the JAVELIN method of Zu et al. (2011), τ JAV = 3.68±0.36 days (Grier et al. 2013a). The transfer function is mainly contributed by the second Gaussian component centered at τ = 3.33 days.
On the top panel of Figure 9, we compare the transfer function derived in this work with those derived from the RLI method (Skielboe et al. 2015) and from the BLR dynamical method (Li et al. 2013). Our result well matches that from the dynamical modeling method but differs from the RLI result, which has a plateau from 0 to about 7 days. The shape of our obtained transfer function also coincides with that from the maximum entropy method (see Figure 1 in Bentz et al. 2010).

SBS 1116+583A
As we have done for for Arp 151, we use the B−band photometric light curve for this object. The range of the time lag in the calculations is set to (0-10) days, and the best number of Gaussian functions is K = 5 according to the AICc. The results are shown on the middle panels of Figure 8. The obtained mean time lag is τ m = 2.50 ± 0.70 days, which is again well consistent with the results based on CCF analysis, τ CCF = 2.38±0.58 days ), and the JAVELIN method, τ JAV = 2.46 ± 0.93 days (Grier et al. 2013a).
The middle panel of Figure 9 compares the derived transfer functions from various methods. The RLI method yields a bimodal transfer function with a response component of around 2 days and an additional component of around 8 days (Skielboe et al. 2015). However, both the transfer functions obtained in this work and those from the dynamical modeling do not show the second component. An inspection of the reconstructed Hβ light curves (see Figure 10 in Skielboe et al. 2015) indicates that our results better match the observations.

Mrk 817
Mrk 817 was monitored for as long as five years by Peterson et al. (1998), but only during the intervening three years were Hβ time lags detected. We use the first-year section of the light curves with the Hβ time lags detected (HJD 24,449,449,211). The CCF analysis by Peterson et al. (1998) yielded a time lag of τ CCF = 20.1 ± 4.4 days. The range of the time lag in our calculations is set to (0-50) days and the best number of Gaussian functions is found to be K = 4. It is apparent from a simple visual inspection that the light curves of the continuum and Hβ line show different long-term trends, with that of the continuum fluxes decreasing more steeply. We therefore include the detrending procedure as described in Section 2.2, in which both the trends of both the continuum and the Hβ line are modeled by a linear polynomial. The results are plotted on the bottom panels of Figure 8. The transfer function peaks at around 17 days and bears large uncertainties due to poor data sampling. The resulting mean time lag also has a fairly large uncertainty τ m = 16.7 ± 8.3 days. Nevertheless, all the variation features in the light curves are well reproduced.
This object has never been analyzed by the RLI method or the maximum entropy method. The bottom panel of Figure 9 only shows a comparison of our result with only that of the dynamical modeling method. The two results are generally consistent, although our result peaks at a slightly shorter time lag.
In our application to the above three RM objects, the transfer functions can be effectively modeled by a single Gaussian. As illustrated in Section 3.2, the time resolution of the obtained transfer function is determined by the data sampling and measurement noises of the observations. In a statistical sense, based on the AICc criteria, the present solutions are the best ones that balance the trade-off between the goodness of fit and the complexity of our approach. Refining the sampling and signal-to-noise ratio can certainly provide more constraints on the fine structures in transfer functions (see Figure 2) if they are not intrinsically a simple single Gaussian. 5. DISCUSSION AND CONCLUSIONS By extending the previous work of Rybicki & Press (1992) and Zu et al. (2011), we describe a non-parametric approach to determine the transfer function in RM. The basic principle that the transfer function can be expressed as a sum of a family of relatively displaced Gaussian functions enables us to largely obviate the need for the presumption of a specified transfer function or BLR model as in previous studies. Thus, this to some extent makes our approach model independent. In addition, following the framework of Zu et al. (2011), the inclusion of the statistical modeling of continuum variations as a DRW process allows us to naturally take into account The black line with shaded areas is the transfer function from this work, the blue line with error bars is from the regularized linear inversion (RLI) method of Skielboe et al. (2015), and the red line is from the BLR dynamical modeling method of Li et al. (2013). Note that there is no RLI solution for Mrk 817. All the transfer functions are normalized on a common scale.
the measurement errors and make use of the information in observation data sets as much as possible. The application of our approach to RM data sets illustrates its fidelity of our approach and its capability to deal with long-term, secular variations in light curves. In particular, our approach is apt for solving the multimodal transfer functions associated with multiple-component structures of BLRs (e.g., Hu et al. 2012).
In addiation to BLRs, our approach can also be applied to RM observations of dusty tori in AGNs (e.g., Suganuma et al. 2006;Koshida et al. 2014).
Since we relax the presumption of a specified transfer function but let its shape be determined by the data, our approach can be regarded as a useful complementary to the BLR dynamical modeling (Pancoast et al. , 2014aLi et al. 2013). Indeed, interpretation of the obtained transfer functions has to invoke physical BLR models. A combination of these two methods may help to remove the concern over the systematic errors in the dynamical modeling. On the other hand, in our approach (compared with the RLI method and the maximum entropy method), we incorporate the statistical description for the light curves of both the continuum and the emission line; therefore, we take full advantage of the information in the data and can estimate uncertainties for the obtained transfer function in a self-consistent manner.
We end with discussions on some potential improvements to be made to our approach.
• Throughout our calculations, we have by default assumed a linear response of the emission line to the continuum so that the covariance functions of the light curves can be analytically expressed by the error function. To include nonlinear response (Li et al. 2013), we can introduce a new parameter γ for the emission line fluxes as so that the new flux series f l,new (t) now linearly responds to the continuum.
• Since our approach is based on the statistical modeling of the temporal variations, it is not straightforward to directly include the velocity information. We can apply our approach to each velocity bin individually (e.g., Skielboe et al. 2015). In this case, any correlation information between velocity bins is neglected; thus, the obtained velocity-resolved transfer function may be noisy along the velocity direction. It merits a future development of a feasible scheme to self-consistently take into account the velocity information.
• We use the DRW process to describe the continuum variations, a process found to be sufficiently adequate for most of current AGN variability data sets (e.g., Kelly et al. 2009;MacLeod et al. 2010;Li et al. 2013;Zu et al. 2013). However, there is evidence that the DRW process is no longer appropriate on a short time scale of days (Mushotzky et al. 2011;Kasliwal et al. 2015;Kozlowski 2016). This impedes the application of our approach in recovering very fine structures (on a timescale of days) in the transfer function. Recently, Kelly et al. (2014) developed a flexible method that statistically describes the variability by a general continuous-time autoregressive moving average (CARMA) process. The DRW process is a special case of the CARMA process. In principle, our approach can also naturally include the CARMA process, but this will make out approach highly computationally expensive. We therefore defer this to a future improvement with the parallelization of our approach on supercomputer clusters.
We thank the two referees (astronomer and statistician) for their insightful suggestions, which improved the manuscript. We are grateful to the staff of the Lijiang station of the Yunnan Observatories for their assistance with the RM data sets. YRL thanks Ying Zu for useful discussions and Andreas Skielboe for providing the data of the Hβ light curves and their derived transfer function for SBS 1116+583A. This research is supported in part by the Strategic Priority Research Program -The Emergence of Cosmological Structures of the Chinese Academy of Sciences, Grant No. XDB09000000, and by the NSFC (Gran No. NSFC-11173023, -11133006, -11233003, -11303026, -11573026, and -U1431228). This work is also supported in part by grant 2016YFA0400700 from the Ministry of Science and Technology of China. This work has made use of data from the Lick AGN Monitoring Project public data release.
We developed for the present approach a software package MICA (Multiple and Inhomogeneous Component Analysis), which is available upon request.

APPENDIX A COVARIANCE FUNCTIONS
In this appendix, we derive analytical expressions for S lc and S ll in Equations (4) and (5). In deriving S lc , for the sake of clarity, we denote where S k lc represents the contribution to the covariance from the kth Gaussian component. It is easy to show where ∆T 1 = ∆t − τ k and erfc(x) is the complementary error function. Similarly, we denote The covariance component S km ll is Introducing new variables u = (τ + τ ′ ) and v = (τ − τ ′ ), we have Note that The above integral can be rewritten as The first integration term is and the second integration term is where ∆T 2 = ∆t − (τ k − τ m ). As a result, we arrive at S km ll (∆t) = πω 2 f m f k σ 2 e ω 2 /τ 2 d e −∆T2/τd erfc − ∆T 2 2ω

APPENDIX B ERROR FUNCTION
The error function is defined as The complementary error function is defined as It is trivial to show that the integral and the intergal e −(at 2 +bt+c) dt = − 1 2 π a e (b 2 −4ac)/4a APPENDIX C COMPARISON WITH JAVELIN The present approach extends the work of Zu et al. (2011) by adopting a more flexible transfer function rather than a specified top-hat transfer function. Therefore, we can obtain more than the characteristic time lag and recover complicated structures in the transfer function. In Figure 10, we run the public package JAVELIN, developed by Zu et al. (2011), on the simulated light curves of the three tests in Section 3.1 and plot the obtained top-hat transfer functions and the distributions of the mean time lag of the top-hat. The recovered transfer functions from our approach and the input transfer functions are also superimposed for comparison. JAVELIN also employs an MCMC method (but with a different algorithm) to determine the best values of the free parameters. Again, we run the Markov chain with 150,000 steps. The time-lag distribution for the first test is bimodal, with the two peaks at 15.53±0.04 and 16.5±0.1 days, respectively, and the overall mean lag of 16.2±0.5 days. Here, the uncertainties are assigned the standard deviations of the distributions. The two peaks correspond to neither the mean lags of the input two top hats (7.5 and 17.5 days) nor the overall mean lag (14.2 days). The JAVELIN time lags for the second and third tests are 14.0±0.2 and 12.5±0.1 days, respectively, in agreement with the mean lags of the input transfer functions (14.0 and 12.5 days). However, for the second test, the time lag does not show the expected double-peak distribution. By contrast, the present approach works fairly well in all these three cases (see Figure 3), in terms of recovering both the shape of the transfer functions and the characteristic time lags. The present approach yields a mean time lag of 14.1±0.5, 14.1±0.9, and 12.4±0.5 days for the three tests, respectively, generally consistent with the input values. Here, the uncertainties are calculated by applying the error propagation formulae to Equation (22). We attribute such discrepancies of JAVELIN to the too simplified top-hat transfer function adopted in the approach, which may encounter difficulties when the transfer function is multimodal and asymmetric.