Regularization methods used in error analysis of solar particle spectra measured on SOHO/EPHIN
A. Kharytonov  E. Böhm  R. F. WimmerSchweingruber
Institute of Experimental and Applied Physics, Kiel University, Leibnizstrasse 11, 24118 Kiel, Germany
Received 25 August 2008 / Accepted 9 December 2008
Abstract
Context. The telescope EPHIN (Electron, Proton, Helium INstrument) on the SOHO (SOlar and Heliospheric Observatory) spacecraft measures the energy deposit of solar particles passing through the detector system. The original energy spectrum of solar particles is obtained by regularization methods from EPHIN measurements. It is important not only to obtain the solution of this inverse problem but also to estimate errors or uncertainties of the solution.
Aims. The focus of this paper is to evaluate the influence of errors or noise in the instrument response function (IRF) and in the measurements when calculating energy spectra in spacebased observations by regularization methods.
Methods. The basis of solar particle spectra calculation is the Fredholm integral equation with the instrument response function as the kernel that is obtained by the Monte Carlo technique in matrix form. The original integral equation reduces to a singular system of linear algebraic equations. The nonnegative solution is obtained by optimization with constraints. For the starting value we use the solution of the algebraic problem that is calculated by regularization methods such as the singular value decomposition (SVD) or the Tikhonov methods. We estimate the local errors from special algebraic and statistical equations that are considered as direct or inverse problems. Inverse problems for the evaluation of errors are solved by regularization methods.
Results. This inverse approach with error analysis is applied to data from the solar particle event observed by SOHO/EPHIN on day 1996/191. We find that the various methods have different strengths and weaknesses in the treatment of statistical and systematic errors.
Key words: Sun: particle emissi on  methods: data analysis  interplanetary medium  methods: numerical
1 Introduction
In many cases the measurements of a physical process is described by the Fredholm integral equation of the first kind
where f is the unknown and the kernel K and the righthand side z are given functions. The usual approach to solving this inverse problem is the discretization and formulation as a numerical algebra problem (Hansen 1998; Neumaier 1998; Tikhonov & Arsenin 1977). We have a system of linear algebraic equations that is often singular or illconditioned
where is a matrix or discrete forward operator with elements a_{ij}, which describe the measurement process, is the unknown vector with components f_{i}, and is the known vector with measured components z_{i}.
When solving this discretized inverse problem with noisy data one must perform the error analysis. For this we consider two approaches.
 1)
 In the statistical approach, we find the error propagation or the covariance matrix from statistical equations. Errors are as a rule in the vector z.
 2)
 In the algebraic approach we find errors from algebraic equations. Errors (called perturbations in this case) can be both in the vector z and in the matrix .
Error propagation for solving inverse problems is considered in Hansen (1998), Neumaier (1998), and van Wijk et al. (2002).
After Wilkinson (1965), the perturbation theory was developed in numerical linear algebra, covering linear systems, the matrix inverse, the least square problem, and singular value problems. One can refer to all those problems as the error estimation. Most of those results are discussed in Wilkinson (1965), Björck (1991), Higham (1996), Lawson & Hanson (1974), Neumaier (1998), Wedin (1973), and they also contain the information about error bounds for the case when the matrix is square or nonsingular. The error estimation for the singular matrix and discrete inverse problems is considered in Neumaier (1998) and Hansen (1998). Lefebvre et al. (2000) present a method of error propagation for matrix inversion, applicable to a non singular matrix A, including errors in the matrix A. This method may be extended to illposed problems, but is not investigated here. The detailed mathematical and physical description of the solar particle spectra computation from the SOHO/EPHIN data have been presented in Böhm et al. (2007), showing the value of the method of regularization, and a comparison of various methods. Missing there was the estimation of uncertainties in the intensities. Thus, in this paper, we focus on the analysis of uncertainties in general and applying various methods to data obtained with the SOHO/EPHIN sensor. Therefore, we give only a brief formulation of this problem in Sect. 3. In Sect. 2 we summarize the telescope EPHIN details, its modeling and experimental data used in the paper. We focus on the error propagation in Sect. 4.1 and consider the error estimation in Sect. 4.3. The statistical approach to the error analysis in perturbed equations is used in Sect. 4.4. The numerical results of the solar energetic particles spectra and errors in solutions obtained by various methods are represented in Sect. 5.
2 Description of EPHIN modeling and measurements
2.1 Instrument and Monte Carlo simulation
The Electron Proton Helium Instrument (EPHIN), part of the Comprehensive Suprathermal and Energetic Particle Analyzer (COSTEP) (MüllerMellin et al. 1995), on ESA's Solar and Heliospheric Observatory (SOHO) measures solar energetic electrons, protons, and alpha particles in the energy range 150 keV to >5 MeV (electrons) and 4 MeV/n to > 53 MeV/n (protons and helium). EPHIN consists of a stack of five silicon detectors, surrounded by an anticoincidence shield, a plastic scintillator and sixth silicon detector to distinguish between stopping and penetrating particles (Fig. 1). Two passivated ionimplanted detectors (A and B) define the large (83) field of view with a geometric factor of 5.1 cm^{2} sr. The lithiumdrifted silicon detectors C, D, and E stop electrons up to 10 MeV and hydrogen and helium nuclei up to 53 MeV/n. We will only discuss electrons and protons. The generalization to inclusion of helium ions is straightforward. For data analysis the energy deposit is pulseheight analyzed. Onboard processing includes classification into E and Pchannels, depending on the energy deposit in each detectors AE. Thus, the Echannels measure predominantly electrons, the Pchannels predominantly protons, but with cross contamination which can be appreciable, depending on the electrontoproton ratio. The pulseheightanalyzed data contain information on the energy deposit of single particles in detectors AE in a total of 256 energy bins, and the depth in the sensor head reached (B, C, D, or E). For this work, we have rebinned this energy information into 30 logarithmically equidistantly spaced energy intervals. The sensitivity of each of these energy bins has been modeled as discussed in Sect. 2.2. Together, pulseheight words and energy bin sensitivities allow us to reconstruct the original particle spectra^{}. This paper addresses the accuracy of the energy spectrum calculations achievable with the regularization technique.
Figure 1: Side view of the EPHINsensor, a telescope consisting of 6 Sisolid state detectors AF (in green) surrounded by an active anticonincidence scintillator G (in red). 10 simulated beam electrons of 4 MeV are shown in red as well as a pair of generated secondary gamma rays in blue. 

Open with DEXTER 
Figure 2: Matrix A of geometric factors in dependence of the energy deposit and the energy of the beam particle , with , , , and in the different quadrants, as derived from the MonteCarlo simulation of the EPHIN instrument. The scale in energy 030 is equivalent to a logarithmic scale 0.1 to 100 MeV, which is repeated from 30 to 60, so displaying the full matrix for electrons and protons. We have applied logarithmic binning in the energies and to account for the powerlaw particle spectra that are expected in particle acceleration processes in the heliosphere. 

Open with DEXTER 
We used GEANT3 (GEANT3 2006) and GEANT4 (GEANT4 collaboration 2006) of the CERN program library to simulate the detector response of EPHIN to ions and electrons. The experimental data were analyzed with the PAW package of the CERN program library (PAW 1993), and information about energy resolution and triggering requirements of EPHIN were added at this stage. Final graphics were produced by ROOT from the CERN program library (ROOT 2000). Figure 2 shows the matrix as obtained from the Monte Carlo simulation, subdivided into four parts. and are nearly diagonal energy deposit beam energy, ( ), with some leakage towards lower energy deposits due to insensitive material and other losses. The crosstalk channels and show much lower response than and and are considerably less diagonal. Figure 3 shows the matrix of uncertainties corresponding to the geometric factors as given in Fig. 2. For small geometric factors the relative errors may reach 100%, as demonstrated in Fig. 4, displaying the relative error in the geometric factor for the matrix elements on the diagonal. Included is the interval in which there is a null row vector in matrix A (no entries in the beam energy for any energy deposit interval). dA was derived from a combination of statistical uncertainty in Monte Carlo simulation and considerations about the uncertainty in the instrument, such as uniformity of dead layers, and energy resolution.
Figure 3: Errors on geometric factors , , , and as derived from the MonteCarlo simulation of the EPHIN instrument. The description is analogous to Fig. 2. Note the different zaxis scales of Figs. 2 and 3. 

Open with DEXTER 
Figure 4: Relative errors of geometric factors for the diagonal elements of matrix A, null row vector in matrix A means (no entries in the beam energy for any energy deposit interval) (see Figs. 2, 3). 

Open with DEXTER 
2.2 Data used in analysis
Figure 5: Time series of Echannel ( upper panel) and Pchannel ( lower panel) count rates during the solar particle event of days 188 to 192 in 1996. The various curves (red, green, black and blue) show count rates in the energy intervals, as recorded by the SOHOEPHINsensor. Day 1996/191, 1617 h, has been chosen for a detailed analysis: a 24h preevent period on day 1996/188 has been taken as solar and nonsolar interplanetary event background measurement. 

Open with DEXTER 
The solar energetic particle event of day 191 in 1996 has been chosen in this analysis, to illustrate again the regularization method in this work, but with putting main effort into uncertainty estimation. We have analyzed other time periods (Böhm et al. 2007; and GòmezHerrero et al. 2006), but will limit ourselves here to the time period defined above in this work.
Figure 5 shows the count rate in that time period for the Echannels,
supposed to be mainly electrons, and the Pchannels, supposed to be mainly protons.
Day 188 has been taken to represent the background for quiet times and has been used
to correct the solar event count rates. Error analysis has been performed for day 191, 1617 h. Figure 6 shows the total pulse height distribution in the E and Pchannel,
which is used as the zvector in Eq. (2). For the uncertainty, ,
we use the rule of quadratic addition of errors (see Bock & Krischer 1998) in the optimization approach. The terms
are the purely statistical errors associated with each data point (omitting uncertainties in the matrix elements of A).
where is the error in the experimental data day 191, 1617 h, in the event, and is the error in the background observation (day 188) (see Fig. 5).
Figure 6: Measured pulse height (energy deposit) distribution on day 191/1996, 1617 h, all Echannels have been combined and all Pchannels have been combined too. The background as measured on day 188/1996 has been subtracted (see Fig. 5). The binning is equidistant in logarithmic scale, 10 bins per decade, corresponding to the energy deposit , see Fig. 2. 

Open with DEXTER 
Note that the validity of this background estimate for our investigated time period hinges on the observation that other background sources, especially Xrayinduced background, are normally concentrated around the onset of solar particle events. To reduce the influence of such background, we have chosen a part of the event, which is late in its sequence, and by using the coincidence technique X and rays only contribute to the count rate by chance coincidences.
3 The basic integral equation and its solution
The basic integral equation (Eq. (4)) connects the intensity of the energy deposit
and the intensity of the original beam energy
of particles.
A particle entering the instrument deposits an energy
in the instrument. Because of incomplete charge collection and insensitive material (such as dead layers), is lower than the original energy of the particle .
Since the energy deposit is a stochastic process, the statistical and experimental accuracy limits the resolution in .
In addition, because of the different ionization cross sections of electrons and protons (and of other ions), not only depends on ,
but also on particle species. Other important factors determining
given
are angles, dead layers, detector arrangements, etc., and need to be taken into account when simulating detector response to a beam of particles with energy
and intensity
.
Thus, the intensity of the energy deposit,
,
is generally given by
which is the Fredholm integral equation (Eq. (1)) that needs to be solved for the unknown quantity . Here, , the kernel of the integral equation, is the instrument sensitivity to measure given a particle of energy . The kernel, , is then given by
where is the efficiency of detecting a particle with energy in that deposits energy in the range , hitting the detector in the area element at position (x,y) within the solid angle element from the direction given by and . and need to be calculated using a MonteCarlo simulation that is described below. We will assume an isotropic angular distribution of the beam and a homogeneous (uniform) distribution on the sensor area. Equation (4) may also be interpreted as an algebraic equation if the kernel is known in discrete intervals, just as Eq. (1) may be interpreted as Eq. (2),
where
and
We solve Eq. (6) by the following regularization methods:
 a)
 The iterated Tikhonov method: the iterated Tikhonov regularization method using the scheme
where is the regularization parameter , is the identity matrix, is the adjoint matrix and r_{k} is the residual
with f_{0} = 0 gives the solution of the Tikhonov method, .  b)
 Singular value decomposition (SVD) method: it is known that the minimum norm solution of the least square problem
is given by the vector
and is called a pseudosolution of the system Eq. (6), the matrix is the MoorePenrose pseudoinverse of .This solution, , is completely congruent with the solution of the Tikhonov method , Eq. (9).
 c)
 Optimization with constraints: the nonnegative solution
is obtained by
optimization with constraints
and with (Eq. (3)) desribing the uncertainty in the numerator, . is used as the starting value. Here we use the nonstandard notation for the deviation instead of the standard symbol , to avoid any confusion with singular value of the matrices. For more details see Böhm et al. (2007).
Then Eq. (6) can be rewritten in the form
(12)  
(13) 
or, in matrix form,
where and are the electron and proton distributions, and z^{E} and z^{P} are the measured count rates in the E and Pchannels.
4 Error analysis
In this section we consider the wellknown and newly developed methods for the estimation of errors in the solution of the discrete inverse problems.
4.1 Error propagation
We will use the general error propagation law that is known from the theory of linear transformations of random variables and in multi variate statistics (Brandt 1999).
For m linear functions
of random variables
with the covariance matrix
or in matrix form
one can obtain the covariance matrix
where is the matrix of derivatives or the Jacobi matrix. Furthermore we assume that the errors are small,
Equation (17) is known as the general error propagation law.
If all x_{i} are independent, i.e. if
is diagonal, the variances of y_{i} are given by the law of error propagation in the wellknown form
4.1.1 Tikhonov method
The error propagation problem with the Tikhonov regularization method is presented in van Wijk et al. (2002).
The regularized solution of Eq. (2) by the Tikhonov method is
where is the regularized pseudoinverse
From Eq. (17) one can obtain the variancecovariance matrix of
Let . If all z_{i} are independent, i.e. if C_{z} is diagonal, the law of error propagation has the form
4.1.2 SVD method
Formally we obtain the variance
from the direct problem, using the minimum norm solution of the least square problem
where the matrix is the pseudoinverse of the matrix . The general error propagation law for this case is
Let , then Eq. (23) in the component form is
where s_{ij} are components of .
For all z_{i} independent the law of error propagation is
or in matrix form
where is a matrix with elements m_{ij}=s_{ij}^{2}, and is the unknown vector with components , and is defined as in Eq. (11).
4.1.3 Minimum variance method
Now consider f as the known vector and Eq. (2) as the direct problem of a calculating z. We define the variance
with components
from the general error propagation law Eq. (17)
where and are variancecovariance matrices.
For statistically independent random variables f and z, matrices ,
are diagonal and the law of error propagation is
Equation (29) in the matrix form is
where is a matrix with elements b_{ij}=a_{ij}^{2}, is the unknown vector with components , is the known vector with components . Thus we have system of linear algebraic equations for finding the unknown vector . For this, one must solve the inverse problem with singular matrix . Algebraic and optimization approaches are used for solving this problem. The algebraic solution can be obtained by the SVD method and the iterated Tikhonov regularization method.
In the optimization approach we must find the minimum of a multivariable function.
4.2 Using MATLAB functions in error analysis
In data analysis the MATLAB function fmincon (Matlab 2006) is used for determining the best estimate
for unknown values f from Eq. (11). Information on the errors of
is obtained from the Hessian matrix
with components
as the elements of the symmetric matrix of second derivatives of the minimization function.
The covariance matrix of
is
where the numerical value k_{QL} = 1, when the function being minimized is a sum of squares (Eq. (11)). The square roots of the diagonal elements are the symmetric errors
The MATLAB function lscov allows us to obtain standard errors and the covariance matrix of the function f.
4.3 Error estimation
We first consider the sensitivity of the solution of Eq. (2) to the perturbations in the righthand side z. If z is changed to
then f is changed to
and Eq. (2) is
Subtracting Eq. (2) from Eq. (34), we have
The solution of Eq. (35) is the perturbation
where is the pseudoinverse of matrix .
Now consider the effect of perturbations in
and z. For the perturbed
solution
we have the equation
where and . Then the solution of Eq. (37) is
where is the pseudoinverse of matrix .
The value
from Eq. (37) or the difference between the perturbed solution
and the unperturbed f can be computed by various means. We can write one solution as
Another expression for is presented in Lawson & Hanson (1974)
The case when is considered in Björck (1996). The perturbed solution satisfies the normal equations
Subtracting , and neglecting secondorder terms, one finds
where .
Wedin & Wikström (2002) use the following formula
where the first term
correspond to perturbations in and the second one, , to perturbations in z (see Eq. (35)). We use Eqs. (36), (42), and (43) in uncertainty calculations described in this subsection.
4.4 Monte Carlo method with random perturbations (Method of numerical ``measurements'')
We consider the perturbations in
and z as a random process. For our case we have
.
The matrices
and
are nonnegative (
and
). Let Q be the number of calculations that one can
consider as the number of ``measurements'' or the sample of a size Q. Then for every calculation k from this set Q the random perturbations matrix
is obtained by
the following method. We use the random matrix
with components
r_{ij}^{k} and
.
Random numbers in this matrix have a normal distribution with mean 0 and standard deviation 1. Then we obtain the matrix
with components
as
where are components of the matrix . Then we form the matrix
with the necessary condition that these new matrices must be nonnegative. We obtain the random perturbations vector by the similar method
where Z^{k} is an array of random numbers with mean 0 and standard deviation 1 that has the same size as z. The vector z^{k} must be nonnegative
If and z for every k (k = 1,...,Q) are changed to and z^{k} then the solution f is changed to .
We solve the algebraic equations
and the optimization problem as the minimum of a constrained nonlinear multivariable function for every k (k = 1,...,Q)
Here z^{k}_{i} are the components of , a^{k}_{ij} are the components of , and are given numbers.
Figure 7: Energy distribution for day 191/1996, 1617 h, solutions and are obtained by applying Eqs. (10), (11). Null row vector in matrix A means: no entries in the beam energy for any energy deposit interval, creating unphysical results, but still producing mathematically correct numbers, small positive or even large negative numbers and uncertainties (see Figs. 2, 3). Points at 10^{2} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER 
Solving these problems we obtain 2Q solutions.
The arithmetic mean of all Q calculations for the function f is
We define the sample standard deviation s(i) for each point i of the solution f as
The error or uncertainty in the mean value in general is the standard error, s_{m}(i), which is defined as
Here, we use this statistical method to get the distribution of values f_{i} with uncertainties in z_{i} and/or a_{ik} that is we use
We can quote the resultant calculations of the solar energy spectra as
where .
This method allows us to consider random perturbations either only in the matrix or only in the vector z or in and z. Besides we can combine the method of random perturbations with the error propagation law. In this case we must use Eq. (24) for every k from Q computations and then evaluate mean values for all three variants of perturbations.
5 Numerical results
Figure 7 shows the energy distribution corresponding to the day 191/1996 event obtained by the regularization methods using algebraic (Eq. (10)) and optimization (Eq. (11)) approaches. These methods produce different and incompatible solutions partly related to the constraints applied in the optimization approach (see Eq. (11)).
Table 1: Summary of error analysis methods.
The algebraic approach produces underflows, small positive values, and even (unreasonably) large negative and positive values. These values are partly due to bins with z_{i}=0 and(or) a_{ik}=0 in a complete column or row (null row vector, see Fig. 7) (see Figs. 24). Underflows in the optimization approach are consistent with f_{i}=0. These energy distributions from Fig. 7 are used in the following as the reference, when displaying relative errors.
Figure 8: Absolute of relative errors in the energy distribution obtained by various methods: the Tikhonov method with , Eq. (21); the SVD method, Eq. (24). Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval). (see Figs. 2, 3) Points at 10^{3} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER 
Figure 9: Absolute of relative errors in the energy distribution obtained by using the minimum variance method, Eq. (30), with using algebraic and optimization approaches. The SVD method, Eq. (24), is placed for the comparison. Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval). (see Figs. 2, 3) Points at 10^{3} display underflows. Points at 10 display overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER 
Figure 10: Absolute of relative errors in the energy distribution obtained by using the MATLAB functions: fmincon, Eq. (33), , and lscov, . Points at 10^{3} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER 
Figure 11: Absolute of relative errors in the energy distribution obtained by using errors in z and ( ) (Eq. (42)), in ( ) (Eq. (43)), and in z ( (Eq. (35)). Points at 10^{3} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER 
Figure 12: Absolute of relative errors in the energy distribution obtained by using the error propagation method (see Sect. 4.1) with random perturbations: in z and ( ), in ( ), and in z ( ), Eq. (48). Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval). (see Figs. 2, 3) Points at 10^{3} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to , see Fig. 2. 

Open with DEXTER 
Figure 13: Absolute of relative errors in the energy distribution obtained by using the Monte Carlo method (see Sect. 4.4) with random perturbations: in z and ( ), in ( ), and in z ( ) for the algebraic approach, Eq. (48). Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval). (see Figs. 2, 3) Points at 10^{4} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to , see Fig. 2. 

Open with DEXTER 
Figure 14: Absolute of relative errors in the energy distribution obtained by using the Monte Carlo method (see Sect. 4.4) with random perturbations: in z and ), in ( ), and in z ( ) for the optimization approach, Eqs. (49), (50). Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval) (see Figs. 2, 3). Points at 10^{4} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER 
Figure 15: Energy spectrum obtained by the Monte Carlo method (see Sect. 4.4), algebraic solution, with uncertainties in z and . errors taken from Eq. (54). The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER 
Figure 16: Energy spectrum obtained by the Monte Carlo method (see Sect. 4.4), optimization solution, with uncertainties in z and . Errors taken from Eq. (54). The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER 
We consider solutions and the error analysis for Eq. (2). Errors in the vector z and in the matrix cause errors in the solution f. Poissonian errors go with the vector z, the experimental data. Errors in the matrix are the result of the Monte Carlo modeling of processes, the modeling of the sensor head, and the statistical accuracy in the sensor response to the particle beam from the Monte Carlo simulation. In all following calculations we define the changes in the solution f caused by statistical errors in experimental data and errors in the matrix represented in Fig. 3. In all results obtained by the optimization approach we use errors only ( in Eq. (11)).
In the following, we compare the various methods for error estimation outlined in Sect. 4 combined with the analysis method explained in Sect. 3. Table 1 gives an overview of the different combinations and methods investigated here and maps them to the figures in this paper.
Figure 8 shows the relative errors in the regularized energy distribution from the error propagation law Eq. (17) obtained by the Tikhonov method and the SVD method. The solution for the Tikhonov method with the regularization parameter is the same as for the SVD method. The errors of the methods coincide where they do not underflow. The errors increase with channel number that is with energy for electrons as well as for protons, and are all below 1 (100%). The SVDmethod produces reasonable errors even in the ``null vector range'', which should not be meaningful.
Figure 9 presents solutions of the inverse problem for error propagation Eq. (30) using algebraic and optimization approaches. For this case the relative errors from the optimization approach are smaller than from the algebraic approach. The errors are on a 10% level, increasing with energy (channel number).
Figure 10 is a plot of relative errors in the energy distribution obtained by using MATLAB functions. The function fmincon allows us to compute the Hessian matrix that is connected with the covariance matrix by Eq. (32). The function lscov allows us to obtain the least square solution in presence of known covariance. We use the following variant of this function (that means optimization with ), here mse is the mean squared error, S is the estimated covariance matrix of f. Errors below channel number 15 for electrons, and number 55 for protons look reliable, staying below 1 (100%). The results with lscov here are not a good measure of the errors ( is used), and should be taken as an example of bad adjustment of the optimization approach with .
Figure 11 shows the relative errors in the energy distribution obtained by using perturbations in the matrix and vector z (Eq. (42)). The main contribution to the errors is produced by errors in the matrix A. It would be possible to improve the calculation of the matrix A, from a purely statistical point of view, but that does not help to avoid systematic errors in modeling the geometry, the materials and in the energy calibration of the sensor head, and in modeling interaction processes in the Monte Carlo program GEANT4. Therefore, we have retained the current uncertainties in A as representative of other systematic unknowns.
Figure 12 illustrates the error propagation law with random perturbations both in the matrix and vector z. The results show that the measurement errors play only a secondary role compared with the uncertainties in the matrix . The relative errors from this figure are the same as the relative errors from Fig. 8. Again errors are obtained in the ``null vector range'' range, here for all variants in calculations.
Figure 13 shows the relative errors obtained by using the Monte Carlo method with random perturbations in the matrix and vector z for the algebraic approach (Q = 100, see Eqs. (51), (55)). Large errors are caused by uncertainties in the matrix A in a wide region.
In Fig. 14 we plot of relative errors obtained by using the Monte Carlo method with random perturbations in the matrix and vector z for the optimization approach (Q = 100, see Eqs. (51), (55)).
The results presented in Figs. 13, 14 were obtained for relative errors and . Here is the largest singular value of , z_{2} is the Euclidean length of a vector z. At low channel numbers the uncertainty in f is dominated by the matrix A, where threshold fluctuations in determination of the geometric factor are important, while at high channel numbers the statistical uncertainties in the experimental date z are dominant.
For the algebraic approach (see Fig. 13) the uncertainties in the matrix play the larger role than in the vector z. For the optimization approach (see Fig. 14) the uncertainties in the vector z play the larger role than in the matrix .
In the case of relative errors and for the optimization approach the uncertainties in the vector z play the insignificant role compared with errors in the matrix . These results are not presented in this paper.
Figures 15 and 16 show our final results. Energy distributions (rate/bin/hour) are transformed to energy spectra ( ) and included are the errors according to the MCperturbation method using and (Eq. (3)). The algebraic approach (Fig. 15) seems extremely unreliable, according to the errors, across the energy range >1 MeV, the optimization approach (Fig. 16) displays large errors at high energies, as expected due to low event rate z and low accuracy in the matrix A (because of low acceptance).
6 Discussion
The uncertainties df_{i} in the spectrum f_{i} depend on the statistical accuracy of the data and on the uncertainties in the matrix of geometrical factors. Qualitatively the statistical uncertainties increase with increasing energy due to the generally steep energy spectrum and reduced count rate. But the detection efficiency is small at low energies, despite the high rate, giving rise to large uncertainties, and the efficiency is low at high energies, because of triggering conditions with anticoincidence, again with small geometrical factors and thus large uncertainties. So we expect large uncertainties from the matrix at low energies, and from the matrix and the data (low statistics) at high energies, especially for electrons.
All distributions of relative uncertainties show this behavior, a decrease at low energies (channel number) and a rise after a minimum, but that does not continue to high energies in the SVD and the Tikhonov approaches, instead staying constant (Figs. 8, 9).
There are energy intervals giving unreasonable (unphysical) intensities and uncertainties, which nevertheless produce reasonable absolute uncertainties, from a mathematical point of view, e.g., in the algebraic approach (Fig. 11), and even in the error propagation method (Fig. 12).
In the optimization approach with constraints (intensity nonnegative) the uncertainties get larger than 100% at high energies, indicating that the error propagation method is not applicable any more. The contribution of statistical uncertainties in the experimental data seems to be less important than the uncertainties in the matrix elements, as can be seen in all methods of error analysis, applying the MC(dz, dA) technique (Figs. 1114). This gives a hint to the fact that the response matrix may show systematic errors, despite producing the matrix to the best of our sensor knowledge. Approximations used in the construction of the geometrical model of the EPHIN sensor, as well as inaccuracies in energy calibration and energy resolution may be the cause of those systematic errors.
The energy spectrum obtained by the algebraic approach shows large errors, when using the MC(dz, dA) technique. These are caused by partly negative values in the intensities f_{i} at higher energies (Fig. 15). This deficiency of the algebraic approach does not appear in the optimization approach (Fig. 16), which gives more realistic values of the uncertainties in the regularization method of data analysis in this experiment.
7 Conclusions
In this paper we have studied the impact of various error sources on the quantity of uncertainties in the solution of the illposed problem Eq. (2). We considered errors in the experimental data (vector z) and the discrete forward operator (matrix A). In the error analysis we use wellknown and newly developed methods, such as the inverse problem for error propagation and the Monte Carlo method with random perturbations in the matrix and vector z. The results for the various methods are presented in this paper.
The algebraic methods, Eq. (10), in general show unphysical negative values in some components f_{i} of the solution f, and in some components of as well. The optimization approach seems to be reliable and stable. The error estimation obtained here correlates with the statistical accuracy in the components z_{i} and the discrete forward operator A.
We find that the various methods have different strengths and weaknesses in the treatment of statistical and systematic errors. Using the MonteCarlo method with random perturbations we managed to include errors in both, vector z as well as in matrix A in our analysis. Based on results shown in Figs. 15 and 16 we claim that this MCmethod provides the best error estimation.
Acknowledgements
We wish to thank the anonymous referee for comments that helped us to improve this paper.
References
 Bevington, P. R. 1969, Data Reduction and Error Analysis for Physical Sciences (New York: McGrawHill) (In the text)
 Björck, A. 1991, BIT, 31, 238 [CrossRef] (In the text)
 Björck, A. 1996, Numerical Methods for Least Squares Problems (Philadelphia: SIAM) (In the text)
 Bock, R. K., & Krischer, W. 1998, The Data Analysis BriefBook (Berlin: SpringerVerlag) (In the text)
 Böhm, E., Kharytonov, A., & WimmerSchweingruber, R. F. 2007, A&A, 473, 673 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Brandt, S. 1999, Data Analysis: Statistical and Computational Methods for Scientists and Engineers (New York: SpringerVerlag) (In the text)
 Branham, R. L. 1990, Scientific Data Analysis: an Introduction to Overdetermined Systems (New York: SpringerVerlag) (In the text)
 GEANT3 2006, CERN Program Library W5013, CERNCN division, see also: http://wwwinfo.cern.ch/geant3 (In the text)
 GEANT4 collaboration 2006, An ObjectOriented Toolkit for Simulation in HEP, CERNLHCC 9844, see also: http://wwwinfo.cern.ch/asd/geant4/geant4.html (In the text)
 GòmezHerrero, R., Klassen, A., Heber, B., et al. 2006, Second Solar Orbiter Workshop, Athens, 1620 October (In the text)
 Hansen, P. C. 1998, RankDeficient and Discrete IllPosed Problems: Numerical Aspects of Linear Inversion (Philadelphia: SIAM)
 Higham, N. J. 1996, Accuracy and Stability of Numerical Algorithms (Philadelphia: SIAM) (In the text)
 Lawson, C. L., & Hanson, R. J. 1974, Solving Least Squares Problems (Eaglewood Cliffs, New Jersey: Prentice Hall, Inc.) (In the text)
 Lefebvre, M., Keeler, R. K., Sobie, R., & White, J. 2000, Nuclear Instruments and Methods in Physics research, Section A, 451, 520 (In the text)
 Lyons, L. 1989, Statistics for Physicists (Cambridge: Cambridge University Press) (In the text)
 Matlab 2006, Optimization Toolbox User's Guide, The MathWorks, Inc., see also: http://www.mathworks.com (In the text)
 MüllerMellin, R., Kunow, H., Fleißner, V., et al. 1995, Sol. Phys., 162, 483 [NASA ADS] [CrossRef] (In the text)
 Neumaier, A. 1998, SIAM Rev., 40, 636 [CrossRef]
 PAW 1993, Paw Analysis Workstation, CERN Program Library Q121, CERNCN division, see also: http://cern.ch/paw (In the text)
 ROOT, Brun, R. 2000, An ObjectOriented Data Analysis Framework, CERNCN division, see also: http://root.cern.ch (In the text)
 Tellinghuisen, J. 2001, J. Phys. Chem. A, 105, 3917 [CrossRef] (In the text)
 Tikhonov, A. N., & Arsenin, V. Y. 1977, Solution of IllPosed Problems (New York: Wiley)
 van Wijk, K., Scales, J. A., Navidi, W., & Tenorio, L. 2002, Geophys. J. Int., 149, 625 [NASA ADS] [CrossRef] (In the text)
 Wedin, P.A. 1973, BIT, 13, 217 [CrossRef] (In the text)
 Wedin, P.A., & Wikström, G. 2002, Technical Report, UMINF:02.16, 1 (In the text)
 Wilkinson, J. H. 1965, The Algebraic Eigenvalue Problem (Oxford: Clarendon Press) (In the text)
Footnotes
 ... spectra^{}
 During times of high fluxes, the pulse height words need to be weighted according to another data type, the socalled histograms to account for the limited telemetry to earth.
All Tables
Table 1: Summary of error analysis methods.
All Figures
Figure 1: Side view of the EPHINsensor, a telescope consisting of 6 Sisolid state detectors AF (in green) surrounded by an active anticonincidence scintillator G (in red). 10 simulated beam electrons of 4 MeV are shown in red as well as a pair of generated secondary gamma rays in blue. 

Open with DEXTER  
In the text 
Figure 2: Matrix A of geometric factors in dependence of the energy deposit and the energy of the beam particle , with , , , and in the different quadrants, as derived from the MonteCarlo simulation of the EPHIN instrument. The scale in energy 030 is equivalent to a logarithmic scale 0.1 to 100 MeV, which is repeated from 30 to 60, so displaying the full matrix for electrons and protons. We have applied logarithmic binning in the energies and to account for the powerlaw particle spectra that are expected in particle acceleration processes in the heliosphere. 

Open with DEXTER  
In the text 
Figure 3: Errors on geometric factors , , , and as derived from the MonteCarlo simulation of the EPHIN instrument. The description is analogous to Fig. 2. Note the different zaxis scales of Figs. 2 and 3. 

Open with DEXTER  
In the text 
Figure 4: Relative errors of geometric factors for the diagonal elements of matrix A, null row vector in matrix A means (no entries in the beam energy for any energy deposit interval) (see Figs. 2, 3). 

Open with DEXTER  
In the text 
Figure 5: Time series of Echannel ( upper panel) and Pchannel ( lower panel) count rates during the solar particle event of days 188 to 192 in 1996. The various curves (red, green, black and blue) show count rates in the energy intervals, as recorded by the SOHOEPHINsensor. Day 1996/191, 1617 h, has been chosen for a detailed analysis: a 24h preevent period on day 1996/188 has been taken as solar and nonsolar interplanetary event background measurement. 

Open with DEXTER  
In the text 
Figure 6: Measured pulse height (energy deposit) distribution on day 191/1996, 1617 h, all Echannels have been combined and all Pchannels have been combined too. The background as measured on day 188/1996 has been subtracted (see Fig. 5). The binning is equidistant in logarithmic scale, 10 bins per decade, corresponding to the energy deposit , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 7: Energy distribution for day 191/1996, 1617 h, solutions and are obtained by applying Eqs. (10), (11). Null row vector in matrix A means: no entries in the beam energy for any energy deposit interval, creating unphysical results, but still producing mathematically correct numbers, small positive or even large negative numbers and uncertainties (see Figs. 2, 3). Points at 10^{2} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 8: Absolute of relative errors in the energy distribution obtained by various methods: the Tikhonov method with , Eq. (21); the SVD method, Eq. (24). Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval). (see Figs. 2, 3) Points at 10^{3} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 9: Absolute of relative errors in the energy distribution obtained by using the minimum variance method, Eq. (30), with using algebraic and optimization approaches. The SVD method, Eq. (24), is placed for the comparison. Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval). (see Figs. 2, 3) Points at 10^{3} display underflows. Points at 10 display overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 10: Absolute of relative errors in the energy distribution obtained by using the MATLAB functions: fmincon, Eq. (33), , and lscov, . Points at 10^{3} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 11: Absolute of relative errors in the energy distribution obtained by using errors in z and ( ) (Eq. (42)), in ( ) (Eq. (43)), and in z ( (Eq. (35)). Points at 10^{3} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 12: Absolute of relative errors in the energy distribution obtained by using the error propagation method (see Sect. 4.1) with random perturbations: in z and ( ), in ( ), and in z ( ), Eq. (48). Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval). (see Figs. 2, 3) Points at 10^{3} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 13: Absolute of relative errors in the energy distribution obtained by using the Monte Carlo method (see Sect. 4.4) with random perturbations: in z and ( ), in ( ), and in z ( ) for the algebraic approach, Eq. (48). Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval). (see Figs. 2, 3) Points at 10^{4} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 14: Absolute of relative errors in the energy distribution obtained by using the Monte Carlo method (see Sect. 4.4) with random perturbations: in z and ), in ( ), and in z ( ) for the optimization approach, Eqs. (49), (50). Null row vector in matrix A means (no entries in the beam energy for any energy deposit interval) (see Figs. 2, 3). Points at 10^{4} display underflows. There are no overflows. The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 15: Energy spectrum obtained by the Monte Carlo method (see Sect. 4.4), algebraic solution, with uncertainties in z and . errors taken from Eq. (54). The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER  
In the text 
Figure 16: Energy spectrum obtained by the Monte Carlo method (see Sect. 4.4), optimization solution, with uncertainties in z and . Errors taken from Eq. (54). The binning is equidistant in logarithmic scale, 10 bins per decade corresponding to the beam energy , see Fig. 2. 

Open with DEXTER  
In the text 
Copyright ESO 2009