Maximum likelihood method for fitting sum of expotentials to experimental dat

Beauregard, R. A.

English

NASA Technical Reports Server

General Disclaimer 
One or more of the Following Statements may affect this Document 
 
 This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 
 
 This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 
 
 This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 
 
 This document is paginated as submitted by the original source. 
 
 Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 
 
 
 
 
 
 
 
Produced by the NASA Center for Aerospace Information (CASI) 
https://ntrs.nasa.gov/search.jsp?R=19670003658 2020-03-16T16:51:57+00:00Z
•nnu
	 +
Hard copy (HCl
MAXIMUM LIKELIHOOD METHOD FOR FITTING
A SUM OF EXPONENTIALS TO EXPERIMENTAL
DATA
By R.A. Beauregard
University of New Hampshire
Department of Physics
Durham, New Hampshire
This work was supported by the National Aeronautics and Space
Administration under Grant #NGR-30-002-018 and Grant #NGR-30-002-021,
and by the University of New Hampshire Computation Center Graduate
Assistantship.
MAXIMUM LIKELIHOOD METHOD FOR PITTING A SUM
OF EXPONENTIALS TO EXPERIMENTAL DATA
R. A. Beauregard
INTRODUCTION
This paper describes a procedure for fitting a sum of
exponentials to a set of experimental data points by the method
of maximum likelihood as described by Qrard [1]'. In particular,
an iterative process is used to maximize the likelihood function.
Part I describes how the likelihood function may be written
when a normal distribution of deviations in the observed data
is assumed. A brief description of the error matrix is given.
Part II describes a modified compute. program (Fortran) designed
for use at the University of New Hampshire using the IBM 1620
and the IBM 360. A sample set of data to be fit by two ex-
ponentials is given along with the solution by the program
described.
[1] The computer program is a modified version of the "Malik"
program by 4rard, Lawrence Radiation Laboratory, University
of California (UCRL-10153)• gee also, "A general Program for
Statistical Analysis Using the Maximum-Likelihood Method (Malik)
Program"; F. arard Nucl. Inat. and MTD's, 34_ 9 242-244 (1965).
PART I
MAXIMUM LIKELIHOOD APPROACH
Ausume n measurements y l .., yi .., yn . are made at times
ti ..ti . . tn in some experiment, and it is desired to fit this
data to a sum of exponentials given by:
f(t)•yt
J•1
aje-bat
 + ak+l
where the a3 and b  are parameters to be determined. The
problem then, is to find values of 
a  
(Jul,••-k+l) amd b^
(j•1, ... K) that gives the best fit of this function to the data.
Further, assume that the deviations of the yi
 from y=f(ti ) are
normally distributed. The deviations are given by:
Ayi • yi - f(ti)
Then the probability of observing a value of yi in the interval
(Y'12'
	
+ 
dyi ) is given by:
1
	
dyi	- (Ayi)2p(yi + dyi) • - .—	 i	 exp	 2-----2	 (1)
i	 i
where (al ) I is the weighting factor to be associated with yi.
The likelihood function is the joint probability of observing
simultaneously the measurements yl••s yi••, y n.•.in the intervals
(Ay149901 009Ayn) and is given by:
dy	 (Ay ) Z
x -=	
exp - i^ 	 (2)L •
161	 '2—s	 °i	 2(°i) _
This assumes all yi measurements are independent.. L in (2) may
be considered as a function of the parameters a  and b, (which
r2-
appear in Ayi ) and we seek values of these parameters such that
L is maximum. To maximise L. it is convenient (and sufficient)
to maximize log L. Hence we define the function F by
n(7j_
 dyi -(Ayi )2	 (3)F • logeL = t	 loge 
	
- Q - -----
	
inl 1	 S	 2(oi)2
It is evident that maximizing F is equivalent to maximizing F +^
given by:
^^-^ ^	 -= 
2
(a)i
The values of the parameters that maximize equation (4) are
the values sought to fit (1) to the observed data. One method
that may be used to maximize ( 4) is described in Part II.
The error matrix associated with the best values of the
parameters may be expressed as fellows. For notational convenience,
we re-label the parameters as follows:
	
41"a1	 42'b1
	
4 3•a2	 44~b2
i This is evidently the "Least Squares" criterion. Note the
assumptions that have been made concerning the deviations in
the measured data to this particular application of the method
of maximum likelihood; that is, the deviations are normally
distributed.
-3-
Let•l o **  2K+1 denote the best values of ¢1••• 12K+1
respectively. We assume that L(41...# 2K+I) is a normal
(multi-gaussian) distribution. Then we may write:
2K+1
L(¢19..9 12K±1 ) ' K exp - 2 E	 A11(¢1 1)(11 ;) (5)3,1x1
where the Matrix 11 Ai1 11 is symmetric, and where the variance-
co-variance of ¢i and ¢1
 is given by:
I I Viol Pit j I' l 0 A11 1 1 -1
To find a convenient expression for the A31 , we proceed as
follows. The log of equation (5) yields:
1 2K+1(2K+1	 ^ ^	 n ^
F•LogL • LogK - - E	 E	 Ai1 ( 1111 _¢i ¢1 -¢^Ii + 4111)
	
2 iml	 301
We now calculate all 2 nd order partials of F.
	
2K+1 2K+1	 n	 ^ TdF	 1	 d [A i1 ( 11 11 — 41 11 -41 11 +41 41 ) J
^—	 E	 E
dTT	 '^ Jul, 	 Jul	 a¢T
2K+1
^ i•1 dCA11(Iili¢iyT
X13 i^^i)J
1 2K+1
Jul
1101
- #1 11 + Ii 43)a CA11 (4141 - Ii11
° IT
r
- A2 a^--[if(111i - X111	 X111 + 11A1)J
T
-4—
where the terms not appearing are zero. Thus:
	
201	 2K+1
dy 
^•'2 i,l [AIT (®i +1 )3 ' 2 1
.1 CAT j (1j-4^	 s (24t -2 )T
	
i.0s	 j¢t
62F	
• 2 Ass 2 A  , n s
^- 2 Ass ) ( 2)s as s
82F	 • — As s	 (6)
Equation (b) given a convenient form of the error matrix. That
Las the variance•co-variance matrix may be calculated with all
2nd order derivatives of the log of the likelihood function.
Therefore:
0! aiaj P1 J I I • j j Ail j j -1
	
^"
-d2F	 -164 1 6J
PART II
COMPUTER PROGRAM
The program described in this paper has been written for
the analysis of experimental data by use of the method of maximum
likelihood. The main program maximizes the likelihood function
by a simple iterative procedure, and also calculates the error
matrix. Three subroutines are used:
a) Matinv to invert a matrix
b) Funct to calculate F;wlogL
e) Input to read and store the observed data
The latter two subroutines must be written by the user only if
he wishes to solve problem other than the one described in
Part I. All programs are written in Fortran Ii.
Iterative Process in Main Program
The method used to maximize P. the log of the likc1Ihood
funetion) is the following. First initial guesses for the para-
meters a,, b  are made and F is calculated. The first parameter
Is increased or decreased by a quantity corresponding to its
assigned step in order to find a larger value of P. If a larger
value of F is found, the new value of this parameter is used and
the same operation is repeated for the next parameter, etc. This
cycle is repeated as long as A increases. If no improvement can
be made wit') the step sizes being used, the steps are all reduced
by a common factor a, and the next iteration is started. The
Iterative process continues until the difference between the
max values of V in the 3th and Q-1) th  iterations respectively
Is less than 0.001.
46
In a given iteration, if a parameter does not yield a
larger value of P after three consecutive attempts, its step
Is reduced by the factor a and the same iteration continues. In
this case, this stop will not be reduced for the next iteration.
Error Matrix in Main ]Program
The error matrix is computed by inversion of the matrix
built with all second derivatives of F at its maximum. The
second derivatives are calculated by differences. In order to
find an adequate a to be used, the width of the likelihood
function L in the direction of each parameter is determined as
being a positive and a negative quantity one must add to each of
the fitted parameters to Lower L by a-1/2 of gilts maximum value
i.e. to lower F by - 2 ) The a's used to compute the second
differences are taken to be a fraction B of the width of L.
B must be read in as initial information. The second derivative
matrix is the pointwise mean of the two secon difference matrices.
The error matrix is the inverse of the second derivative matrix.
During the calculation of the second differences, it is
possible that a larger value of A is obtained, if the iteration
had stopped at a saddle point. In this case the iterative process
Is re-started with the initial steps.
Input Subroutine
In general, the input subroutine is used ones to read the
observed data points yi , ti , vi . These values are used exclusively
by the Funct subroutine.
-7-
Funct Subroutine
In general, the Funct subroutine Is used repeatedly to
calculate A • log L in terms of the parameters and the observed
data points.
The Funct'subroutine that was written for the problem in
Part I calculates:
nt2 ^•1
x	 2
yj'
	
sie 
_bit,
o^
Input Data
The data cards read by the main program are the tollo#'kng:
1) One card with the number of parameters to be fitted.
Format (20 x, 110) This number must not exceed 10.
2) At most ten data cards, one for each parameter (initial
guess) and the corresponding assigned step.
Format (10 x, 2FI0.5)
3) One card containing the factor a used to reduce the
steps at the end of each iteration, the total number of iterations
allowed, and the factor B used to calculate second differences.
Format (20 x, F10.5, I10 9 F10.5)
The main program then calls the Input subroutine. The data
cards read by the Input subroutine that was written for the
Problem in Part I are the following.
48
4) One card containing the number of data points (maximums
In 40). Format (I5).
5) At most 40 cards, each containing the values ti' Yip
of Format (F10.5, 2F10.Z)
SAMPLE SET OF DATA FITTED TO THE SUM OF TWO EXPONENTIALS
Yi•a1e-blti + a2e-b2ti
Yi ti of
26268 .0204 271.94
14636
.1105 168.73
11200 .1561 145.74
8942
.2008 109.49
6465 .2464 96.28
5308 .2913 74.81
4211
.33362 70.19
3171 .3814 62.61
2543 .4264 59,05
2166 .4714 56.37
1718
.5173 51.71
1561 .5622 41.83
1261
.6069 39.91
1172 .6519 38.87
952 .6971 37.40
903 .7425 33.35
74 7 .7874 32.38
793 .8328 32.85
618
.8778 31.71
6 71 .9321 32.99
492 .9682 32.09
607 1.0134 34.65
579 1.0584 37.33
al . 28320 ± 2.62
b i a 6.519 ± 0.042
a2 = 1245 `# 101
b2 w 0.86 ± O.G$


Maximum likelihood method for fitting a sum of exponentials to experimental data

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server