On error propagation in the resolution correction

bUCLEAR

INSTRUMENTS

AND

METHODS

39 (~966)

271-277;

i=) N O R T I I - H O L L A N I )

ON E R R O R P R O P A G A T I O N I N T H E R E S O L U T I O N

PUBI. ISHING

CO.

CORRECTION*

G. QUITTNER

Osterreichische Studiengesellschaft fiir Atomenergie, Ges.m.b.H., Institute of Physics, Reactor Center Seibersdorf, Austria Received 5 August 1965 If experimental functions subject to random fluctuations are replaced by appropriate smoothed versions, the appearance of spurious oscillations in the results of resolution correction can be avoided. The truncation of series expansions of the experimental functions with respect to a suitable system of orthogonal functions is a smoothing operation of this kind. Three of such smoothing operations, for the input, for the resolution- and for the output-function respectively, are introduced and some comments on their use are given.

An analysis, by means of Fourier transforms, of a simplified resolution equation is used to give insight into the respective roles of the frequency spectrum of the true function, of the uncertainties of the input function and of the resolution function, and of the resolution width in determining the amount of true details in the results of measurement. Implications for the proper number of equations to be used in the solution of the resolution equation and for the choice of the resolution width in an experimental arrangement are considered.

1. Introduction

various factors which p r o d u c e or f a v o u r oscillations in the results and which at the same time d e t e r m i n e the a m o u n t o f true structure therein. The r e m a i n i n g sections give the implications o f sect. 2 and 3 for the choice o f the n u m b e r n in the system (2) and o f the resolution widths and c o m p a r e the results o f the present article with those o f other authors.

We consider the p r o b l e m o f restoring the " t r u e " f u n c t i o n f ( y ) from the experimentally given " d i s t o r t e d " function g(x), the distorting action o f the m e a s u r i n g instrument being described by a kernel k ( x , y ) in the equation g(x) =

f

k ( x , y ) f ( y ) dy.

(1)

a

In numerical calculations, this integral e q u a t i o n is often replaced by a system o f linear e q u a t i o n s

2. The smoothing operations used for the control of error propagation. Cut-off numbers m We first define the o p e r a t i o n s :

Operation A gi = ~

kii.]),

i = 1. . . . .

n

(2)

Expand the given input function g(x) into a series

j=l

In the rather extended literature on this subject [see e.g. ~) for a reference list] some difficulties have been r e p o r t e d in a p p l y i n g the respective methods, the ~rouble being generally the a p p e a r a n c e o f spurious oscillations in the results. S o m e a u t h o r s suggest, as a remedy, the use o f s m o o t h e d functions2). Others 3) stress the i m p o r t a n c e o f c h o o s i n g the correct n u m b e r of e q u a t i o n s in the system (2). A l t h o u g h it is c o m m o n l y felt that these difficulties and p r o b l e m s have to d o with the p r o p a g a t i o n o f statistical fluctuations o f the input d a t a into the result, there seems to be no systematic t r e a t m e n t o f this p r o b l e m giving a satisfactory insight into the relevant factors and positive suggestions how to o v e r c o m e the difficulties in actual cases. In sect. 2 we introduce the o p e r a t i o n s for the c o n t r o l o f e r r o r p r o p a g a t i o n in the numerical execution o f the resolution correction and give some c o m m e n t s on their use. The subject o f sect. 2 is s u p p l e m e n t e d and explained by sect. 3, where we use a s o m e w h a t specialized form o f the resolution function to give insight into the * The research reported in this document has been sponsored in part by the U.S. Government. 271

g(x) = ~ G,q~,(x) i=1

(3)

with respect to a system {~0i(x)} o f o r t h o g o n a l functions o f increasing complexity, allowing successively closer a p p r o x i m a t i o n s o f g(x) in the weighted-least-squares sense. [The quantity to be minimized by the choice o f the a p p r o x i m a t i n g function o f o r d e r p is

,VI, =

akq~k(X,) -- g(xt~) /1=1

k

a2(X¢) , P < S (4)

1

where xp are values o f the i n d e p e n d e n t variable at which measured values g(xp) are available; a2(xtl) are the variances o f these measured values, s is the n u m b e r o f measuring points. Also " o r t h o g o n a l " means " o r t h o g o n a l with weights I/a2(xt~)":

'Pi(-':,) (I/a2(xa))q~j(xa) = 6i~.

(4a)

1~1

With this definition o f the q~i the e x p a n s i o n coefficients ak are statistically independent]. This m a y be the

272

(;. QUITTNER

system of trigonometric functions* or the system o f orthogonal polynomials ~) or some other suitable complete system of mutually orthogonal functions. F r o m the given uncertainty Ag(x), calculate the uncertainties AGi. C o m p a r e each AG~ with the corresponding G~ and omit from series (3) all terms with

I

input data

I

f

I

I ag(x)

AGI >=.~ I O~, ,

(5)

f [

where ~ is a number of the order of 1. The number of remaining terms in series (3) is called the cut-(~'number mrs, the index G pointing to function g(x). The output of operation A is the smoothed version o f g ( x )

[ I

f

J

m G

0(.,) = Z ajq, (x)

(6)

L

j=t

J

and the cut-off-number m~.

Operation B

f (x)

f(x)

Perform operation A o f smoothing and determining cut-off numbers on the resolution function k(x,y). Because ofk(x,y) being two-dimensional, the procedure is more complicated, but not different in principle. The output is a smoothed version E(x,y) o f k(x,y) and the two cut-off numbers mKx and mK~..

• f (x)

Operation C Apply operation A to the o u t p u t f ( y ) of the resolution correction. The output is the smoothed version f ( y ) of fO') and the cut-off number m~.

Comments on operations A, B, C The meaning o f the smoothing operations A, B, C and their relevance for the problem of suppression of oscillations in the results of resolution correction can be seen much better after the discussions o f sect. 3. Some remarks however, seem to be useful for understanding already at this point. The series expansions like (3) are necessary for the separation of the experimentally given data into a part which is "'accepted" and another one which is "rejccted". The rule (5) divides the sum in (3) into terms which are " k n o w n " and others which are " u n k n o w n " , i.e., obscured by statistical fluctuations. It will be seen in sect. 3 that this very " u n k n o w n " terms are responsible for the appearance of spurious oscillations in the results, or to put it in other words, the oscillations are the result o f treating " u n k n o w n " terms as if they were " g o o d " data. The choice of ct is to certain extent arbitrary, depending on what one is willing to accept as "certain", a smaller c~ meaning a more conservative choice. In many cases the series {J Gil}, {1Fil} .... are mono*

Linear combinations of the usual trigonometric functions must be used t o ensure the "'orthogonality with weights".

output

Fig. 1. Flov,' diagram illustrating the place of the operations A, B, C within the unfolding procedure. tonously decreasing in i, so that all terms up to a certain i are retained and all with a higher i arc disc a r d e d - t h e r e f o r e the name "cut-off numbers" for 111(;, m K x , mK), , mr...

The flow diagram in fig. 1 shows the place of operations A, B, C within the process of data reduction. The procedure o f "global" smoothing as described above cannot be replaced by simpler smoothing operations like avering over ncighbours ("local" smoothing) because the latter do not provide a control of the degree o f smoothing which is an essential feature o f the present approach and do not give the cut-off numbers m~, m t..... which are essential for the proper numerical solution o f the resolution equation. 3. Factors determining the error of the results of the resolution correction In the discussions o f this section we neglect the variation of the shape ofk(x,y) over the range o f x and take some mean of the different l,(xi,y)'s, for tixed x i, as our one-dimensional resolution function h(x-y). Eq. (1) is then rephtced by

273

ON E R R O R P R O P A G A T I O N IN THE R E S O L U T I O N C O R R I C T I O N

g(x) =

f

b

h ( x - y)f(y) dy.

(7)

tt

[]'his may be a poor approximation for accurate numerical calculations but serves as a model for obtaining relationships which are qualitatively valid for the real, two-dimensional k(x,y)]. Also for convenience we use Fourier integral transforms instead of a discrete series expansion of eq. (3). Now, if G(p), H(p), F(p) are the Fourier transforms ofg(x), h ( x - y ) , f ( x ) of eq. (7), then by virtue of the "Faltung" theorem [ref. 5) p. 24]

G(p) = H(p)*F(p).

(8)

p is the variable ("frequency") in the space of the Fourier transforms. G(p), H(p), F(p)are thecontinuous counterparts of the series of expansion coefficients G,, H;, F, of sect. 2. Now we show, by a series of schematic diagrams in the p-space, that and how the four factors: 1. the uncertainty Ag(x) of the input experimental function, 2. the uncertainty Ah(x-y) of the resolution function, 3. the width of the resolution function, 4. the spectrum F(p) of the true functions determine the amount of details in the result.To seethe influence of the respective factors some of the cases listed below have to be compared with appropriate reference cases which are idealized in some respect. As a special case we can state the conditions for neglecting the influence of factors, in particular t\3r factor 2. and for factor 3. In the drawings, figs. 2-6, we use some simplifying

t.O

tO

0.5

05

0

p

assumptions: a continuously falling G(p); ct = 1 in the decision rule (5); a Gaussian h(x-y); a constant AG(p) and A H(p). Case I (fig. 2). Zero uncertainty of the measured data, perfect resolution, Ag(x) = O, Ah(x- y) = O, h ( x - y) = ~5(x- y). This is a trivial case, there is no resolution correction at all,

r(p) = G(p). Case 2 (fig. 3). Finite accuracy of measured data, perfect resolution. The role of Ag(x) (ref. case l). Rule A of sect. 2 requires that all frequencies higher than the cut-off frequency P should be omitted. The same cut-off frequency P applies also for the result F(p). No resolution correction. Case 3 (fig. 4) Zero uncertainty of the measured function, h(x-y) of finite width but perfectly known. The width of H(p) is inversely proportional to that of the original h(x-y) [ref. 5) p. 523]. In the transition F(p)~ G(p) the amplitude of the high frequency terms is reduced as compared to the low frequency terms. This is the "smearing out" action of a measuring instrument with a finite resolution width. The transition G(p)~ F(p) corresponding to the resolution correction has the reverse effect of magnifying high-frequency terms. In this (hypothetical) case however the true function f(x) can be completely reconstructed from the knowledge of g(x). Case 4 (fig. 5) General conditions, except that the resolution function is perfectly known Ah(x-y) = O.

H ( p ) . ¢onst

1.0 ¸

0.5

0

p

p

Fig. 2. The Fourier transforms G(p), H ( p ) , F ( p ) of the experimental function g(x), the resolution function h ( x - - ) 9 and the true function f O ' ) . Perfect resolution, zero uncertainty of g(x).

274

~;. QUITTNER

1.0

to

05'

05

F!(p) .

const

tO

0.5

(P)

dl 6(p )

0

.

p

o

p

0

P

p

Fig. 3. Same as fig. 2 but finite accuracy . l g ( x ) of the measured data. P is the cut-off frequency.

The role of the width o/" the resolution ./imction (ref. case 2). The width of the resolution function may influence the accuracy of the result in the actual cases where Ag(x) is finite. This is because G(p) is steeper then in case 2 and therefore, other things being equal, the cutoff-rule (5) applies at a lower frequency Q than in case2. The case 4 illustrates the mechanism of error magnification which is of central importance for the problem of the present article: in the process of resolution correction G(p)--*F(p), not only the high-frequency terms but also their errors are magnified (see the

t.O

t.OI

0.5

0.5

P

AF(p)-curve in fig. 5, right section) and appear as oscillations of large amplitude in the result, the amplitude of the low frequency terms remaining relatively small. Fig. 5 also illustrates why operations A, B, C remove the oscillations from the results: they eliminate just these higher frequency terms ofg(x) which are of small amplitude as compared to their own errors and which when restored to their original magnitude bring about large spurious oscillations in the result. The influence of the finite resolution width is howcvcr negligible as long as Q = P and this is the case if

tOt

05

P

Fig. 4. Functions G(p), II(p), F(p) u h e n the resolution function h ( x - - y ) is a G a u s s i a n of finite ~idth, g(x) and h ( x - y) are perfectly k n o ~ n . This is the hypothetical case o f perfect resolution correction.

ON ERROR P R O P A G A T I O N

IN THE RL:SOLUTION

f.o

I.O

275

CORRECTION

1.0"

i /

i! l/ /'~F(p) 0.5

0.5

/

.5

//

•

O , C_p__~~_ o

/

_ p

o

p

p

Fig. 5. Same as fig. 4 but finite accuracy/_lg(x) of the experimental data. Because of the perfectly known Q is the same for and

G(p) F(p).

h(x--y),the cut-off

frequecny

H(p) is not substantially different from a constant within the range OP of interest. The resolution correction can then be omitted despite the finite resolution width. Case 5 (fig. 6)

of G(p) and H(p) respectively (see also appendix 2). We summarize the conditions for possible simplifications of the resolution correction, in terms of cutoff points in the p-space:

The role o[ the uncertainty Ah(x- y) of the resohttion Jimction (ref. case 4) The error of the quotient F(p)= G(p)/H(p) can be

CONDITIO~ CH.

calculated by standard techniques. The application of the cut-off-rule (5) to F(p) may result in a still lower cut-off point R than in case 4. Whether this is the case can be seen quickly from the relative magnitudes of the cut-off frequencies Q and S

If Q ~ P, the resolution correction can be omitted.

CONmTION CDH. If R ~ Q, the uncertainty of the resolution function can be neglected.

The role of the spectrum F(p) of the true functions. The discussions above show that the effect of all other

I 1.0

1.0

l.OI /

~H(p) 0.5"

0.5

/~F(p) /I / /

•_'?~PJ "---____~ 0

Q

p

S

0

General case; finite resolution, unperfectly known

p

0

R

Fig. 6. The cut-off frequencies Q, S, R for G(p),

g(x)and h(x--y).

p

H(p),F(p)are different.

276

(;. Q U I T T N E R

factors 1, 2 and 3 is larger or smaller depending on whether F(p) is fiat [i.e. a large amount of details in .f(y)] or steep [smooth./'(y)'s]. This is true especially for the conditions CH and C D H . 4. Implications for the choice of number of equations in the numerical solution of the resolution equation In this paper w.e are not concerned with a particular method of solving the resolution eq. (I). The circles of operation A, B, C in the flow diagram of fig. 1 are well separated from the box "'resolution correction" and operations A, B, C only modify the functions which are to be fed into, or are the result o f the resolution correction. A class of procedures tbr numerical solution of eq. (1) use as their starting point the system (2) of linear equations:

gi = L kijJ'J, j=l

i = 1. . . . . n

approximating eq. ( 1). The number n of equations must be chosen correctly: too small a number would destroy useful information, too large a number would generate instabilities. This number n is closely' interrelated with our operations A, B, C. We remember that A, B, C give three different cutoff numbers m 6, m~,~ and m~,?.. In practice m6 will often be the smallest a m o n g them. Then it would be natural to take m6 as the number n in the equations, m,; values of ~(.ve) taken at arbitrary m , points can be used as an input to the resolution correction calculation. However, the function ~(x,y) contains some useless "'harmonics" and in order to avoid the appearance o f them in the result one should then smooth ~(.v,)') to a /<(x,)') which contains the same number of ~p~(.v)harmonics as O(x). Another procedure would be to take (mK~,, ml,-~.)values of ~(x,y) and correspondingly ;n~. values of.q(x). The superfluous harmonics coukt afterwards be removed by operation C. The smaller of rib,.,, and m~:~., is the maximum number of equations permitted by an unperfectly known resolution function (appendix 1). 5. Implications for the choice of the ~idth of the resolution function The problem has been stated in a previous paper [ref. ~), sect. 4]. The condition given there corresponds broadly to condition CH. i.e. to the case where the resolution correction can be neglected. Using the discussions of sect. 3 we are now in a

position to see the consequence of the choice of a larger or smaller resolution width in tile more general case. The consequences o f tile choice of a broader resolution function are three-fl)ld: I. The statistical uncertainty of the uncorrected data, Ag(x), becomes smaller. 2. The spectrum G(p) becomes steeper. (The effects of 1. and 2. on the cut-off frequency Q in fig, 5 are opposite). 3. As smaller high-frequency terms have to be restored to their original magnitude the results depend to a larger extent on the accuracy and reliability of the resolution correction. Therefore the shape of" the resolution function must be known with greater accuracy and this accuracy may become the factor limiting the final accuracy. The choice of a particular resolution width is advantageous only if, for the particular measurement, the final cut-off frequency R is larger than with another resolution width. N o w there are cases where an accurate determinalion of k(x,),) is difficult or even impossible. One is then forced to choose a resolution width small enough to meet condition C D H or even CH, at the expense of the statical accuracy of g(x). In other cases where accurate determinations of k(x,)') are possible, much broader resolution functions can be tolerated. 6. Some connections with the work and results of other authors The smoothing of experimentally given functions by representing them by a suitably truncated series expansion is described e.g. by Klepikov and Sokolov4). Their treatment corresponds to our operation A but they are not concerned with the resolution correction. There is a nt.mber of papers reporting the appearance of oscillations in the results, e.g. 5, u) which, however, contain only a few, mostly qualitative and isolated rcmarks on their causes. In sect. 4 we were able to list four relevant factors and give a discussion of their relative importance in producing or favouring spurious oscillations. Moreover, operations A, B, C are means for suppressing the oscillation without an undue loss of useful information. Several authors, e.g. x) state that the use of smoothed functions would help aw)iding the difficulties but do not give definite prescriptions for generating appropriately smoothed functions. The functions .~(x),/<(x,y)./(r) of sect. 2 and k ( x . y ) of sect. 4 are smoothed to a degree which is ,just adapted to the gixen conditions. The paper most similar in the approach to ours seems to be that of W o r t m a n and Cramer ,It."). They

ON E R R O R

PROPAGATION

IN T H E R E S O L U T I O N

use a n-degree polynomial fit to g(x) to calculate the resolution correction as a difference and apply this difference to the u n s m o o t h e d data. n is stated as an input d a t u m to the c o m p u t e r p r o g r a m but unfortunately no indications on the choice of n are given. Our discussions show that as long as the resolution correction is small, this n may be appreciably smaller than m~ without too a large loss in final accuracy. The choice of n < rn~ simply means that not all of the information contained in the experimental data is used for the resolution correction. The use of a small, constant n instead of our variable rna may be a good conservative practice. Instead of separating the smoothing operations A, B, C from the main resolution correction as in our fig. 1, one may try to incorporate the smoothing operations in the main resolution correction calculation. Especially the iterative methods ~) seem to be suitable for this purpose because they actually do the successive approximation of functions just as our operations A, B, C. Terminating the a p p r o x i m a t i o n at a suitable point TM) would be roughly equivalent to performing operation C alone (fig. 1). When sufficiently developed, such a method may offer the advantage of doing the same work as the scheme of tig. 1 with less computational and p r o g r a m m i n g expense. The advantage of the scheme given in fig. 1 is the high degree of control at each stage of computation. In a similar way other simplified procedures may bc t~sed under certain conditions. The decision whether these conditions are met can be done by comparison with the general scheme of tig. 1 by help of discussions similar to that of sect. 3. The author is indebted to Dr. H. Paul for helpful discussions during the completion of the manuscript and to P. Kernthaler and G. Herzog for calculation and drawing of the figures.

Appendix 1 MAXIMUM NUMBER OF EQUATIONS VCHICH CAN BE USED WITII A RESOLUTION FUNCTION

k(x,y)

OF FINITE

ACCURACY A k

Fox and G o o d w i n 3) mention the following theorem: If k(x,y) is a polynomial in x of degree p, then the (n,n) matrix of sample values of k(x,y) is singular for n_>-p+2. The author is indebted to E. Balcar for pointing to a simple p r o o f of this theorem• By hypothcsis k u = k(xi, )'j) = ~_.,P=o a~(y~)x~ with n _>__p + 2. The matrix /,ij can then bc represented as a matrix product

(kll...kl~

277

CORRECTION

/ = (9)

k.~ . . k . .

=

/

•

,,ao(y~)

• a,O',,)o ..o

(i) XI

",-'

•Xn

.x~~-'

N o w in the first matrix on the right side, there is (by hypothcsis) at least one column of vanishing elements and therefore the corresponding determinant is zero. So the determinant of the matrix k u is also zero and the system (2) of equations is not soluble. N o w if we apply to a given k(x,y) operation B to avoid the propagation of Ak into the result and use polynomials as the functions {¢~}, then k(x,y) is represented by a polynomial of degree m ~ x - 1 in x or by a polynomial of degrce rnKy-1 in y. By applying the above theorem (which is of course also valid for y), the smaller of m~x and m~.is the max. n u m b e r of equations.

Appendix 2 RELATIONSHIP BETWI!I:N THE CUT-OFF FREQUENCIES Q ,

R,S

In the case where Q ~ S and therefore the contributions of Ag and Ah to AJ'are comparable, the final cutoff frequency R is lower than either Q or S. We can demonstrate this fact by simple calculation in the special case of a small ~ (in the order of 0.1) because then the standard error propagation formula

(dr~F) ~ = (AG/G) 2 +(,~H/H) ~

(lO)

applies. We evaluate fl)rmula (10) at the frequency Q which by hypothesis is also the frequency S. At this frequcncy AG/G = AH/ H = :~ and by (10), AF/F = ~,"2. The point R, where by definition A F / F = :c is valid, is therefore at a somewhat lower frequency. In the general case Q ¢: S, R will virtually coincide with the lower of Q and S, i.e., usually with Q.

References 1) R. Gold, TID-18304, ANL (1963), 2) G. Kreisel, Proc. Roy. Soc., London, A 197 (1949) 160. 3) L. Fox and E, T. Goodwin, Phil. Trans. Roy. Sot., London A 245 (1952) 50. 4) N. P. Klepikov, S. N. Sokolov, Analysis and Planningof Experiments by the Methodof Maximum Likelihood(Pergamon, 1961). 51 !. N. Sneddon, Fourier Transforms (McGraw Hill, 1951). 0) R. E. Rand, Nucl. Instr. and Meth. 17 (1962) 65. 7) G. Quinner, Nucl. Instr. and Meth. 31 (1964) 61. 8) W. R. Dixon and J. H. Aitken, Can. J. Phys. 36 (1958) 1962. 9) Wortman, Cramer, Jr., Nucl. Instr. and Meth. 26 (1964) 251. 10) Sharsgard, Johns, Green, Rad. Res. 14 (1961) 261.