### Create a StudySoup account

#### Be part of our community, it's free to join!

Already have a StudySoup account? Login here

# Intro To Comput In Phys PHYS 6720

The U

GPA 3.5

### View Full Document

## 25

## 0

## Popular in Course

## Popular in Physics 2

This 34 page Class Notes was uploaded by Ericka Ryan I on Monday October 26, 2015. The Class Notes belongs to PHYS 6720 at University of Utah taught by Carleton Detar in Fall. Since its upload, it has received 25 views. For similar materials see /class/230021/phys-6720-university-of-utah in Physics 2 at University of Utah.

## Reviews for Intro To Comput In Phys

### What is Karma?

#### Karma is the currency of StudySoup.

#### You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 10/26/15

Physics 6720 Cubic Spline Interpolation October 30 2002 Splines and particularly cubic splines are very popular models for interpo lation Historically a spline77 was a common drafting tool a exible rod that was used to help draw smooth curves connecting widely spaced points The cubic spline curve accomplishes the same result for an interpolation problem In these notes we describe cubic splines and discuss their use in interpolation and curve tting Many texts give the algorithms for determining the spline coef cients so we do not do that here 1 Cubic Spline Interpolation We start from a table of points for 239 01n for the function y That makes 71 1 points and n intervals between them The cubic spline interpolation is a piecewise continuous curve passing through each of the values in the table There is a separate cubic polynomial for each interval each with its own coef cients llz 7 03 biz 7 02 cix 7 d for z 6 sh 11 together these polynomial segments are denoted Sz the spline Since there are n intervals and four coef cients for each we require a total of 471 parameters to de ne the spline We need to nd 471 independent conditions to x them We get two conditions for each interval from the requirement that the cubic polynomial match the values of the table at both ends of the interval ya Silt i1gt 1141 Notice that these conditions result in a piecewise continuous function We still need 271 more conditions Since we would like to make the interpo lation as smooth as possible we require that the rst and second derivatives also be continuous SAW 5119507 SL1lt gt 5151950 These conditions apply for 239 12 n 71 resulting in 2n 71 constraints So we need two more conditions to completely x the spline There are some standard choices left to the user 1 natural77 Sam 0 sgqm 0 1 2 clamped Sam f x0 S klwn f xn Other choices are possible if the function is periodic Which is best depends on the application With 471 coef cients and 471 linear conditions it is straightforward to work out the equations that determine them The conditions can be reduced easily to a tridiagonal system with the coef cients 0 as the unknowns Once solved the remaining coef cients are easily determined There are many good codes in the public domain Cubic splines are popular because they are easy to implement and produce a curve that appears to be seamless As we have seen a straight polynomial interpolation of evenly spaced data tends to build in distortions near the edges of the table Cubic splines avoid this problem but they are only piecewise continuous meaning that a suf ciently high derivative third is discontinous So if the application is sensitive to the smoothness of derivatives higher than second cubic splines may not be the best choice 2 Cubic Spline Smoothing When interpolating we start from reasonably exact tabulated values and re quire that the interpolating function pass exactly through the values ln curve tting we generally start with a table of experimental data in which the y val ues are imperfectly known so come with a con dence interval lmy 0 and we require that the tting function come only reasonably close to the data The table of experimental data consists of three columns ziyi 0 with the last value specifying the con dence range commonly one standard deviation ln cubic spline tting we relax the requirement that the spline pass exactly through the points and demand only that the spline and its rst and second derivatives be continuous SHW 549507 5571 519507 51195139 Ski95139 So continuity still leaves 71 3 degrees of freedom They are determined by balancing two opposing criteria 1 The spline must come reasonably close to the data 2 The spline must have low curvature We use a chi square measure to quantify how close the spline comes to the data n 2 X2 Z 2 1 i0 z The numerator of each term measures the discrepancy between the spline and the data point The denominator provides a weight The higher the uncertainty 0 the farther we are allowed to miss the data point If the spline is a good representation ofthe data so deviations are only caused by statistical uctuations the average value of each term in X2 is expected to be about one For reasonably large n the value of X2 is then about 71 i To quantify the curvature we integrate the square of the second derivative g1v1ng lS zl2dz These constraints are contradictory Notice that to make the curvature zero its smallest possible value the spline would have to have zero curvature 7 Le a straight line But that would probably give us a high value of X2 On the other hand we can make X2 equal to zero by having the spline interpolate the points exactly but that would give a high curvature Putting these two constraints together we require that the cubic spline minimize W pXZ lS zl2dz The constant p determines the tradeoff between the two contradictory criteria Putting p 0 removes the constraint on X2 altogether and allows the spline to become a straight line Putting p very large puts high emphasis on minimizing X2 forcing the spline to interpolate without regard to curvature Minimizing subject to the continuity requirements leads again to a tridi agonal system that is easily solved to give the coef cients of the cubics Making p too small allows a high value of X2 Making it too large forces an unreasonably small value of X2 Usually we have to experiment with the choice The best value results in X2 in the range 71 i We should consider smoothing77 to be a poor substitute for tting data to a decent model function that has a fundamental basis When we t experimental data to a model function we are actually testing our understanding of nature When we smooth in order to t to an arbitrary cubic spline we are merely parameterizing an observation but not learning something more fundamental 39sanm 0113 aq3 30 uo31gt30dxa3u u1gt m paxgnbax s 11 1371 0d 037 cared 868998 081 31 380113 31gt am ap 81gtq VOJ 123m mind aq3 31gtq3 ms 3 581m 83 3g 2 3daaxa 3u1gt3suo3 1gt sg awq 1133133 KJGAG 133ng 058 1581 lx Oxo 1 0 39 A x 3 058 2 158 2 d ua33gm aq U123 satqm 0313 383 aq3 8U 31gtggtd333ug Impng 31213ng 1gt uaql 3911quot 39 quot0 g 303 Wi m Sq 1123 aq3 ug sanpm aq3 Sunmmp Sq ultgt 31gtodm3ug 31213ng aq3 ampwua S331 uogqmodJequI a um eq I 39SZZZ39I Emmi sanm 0313 aq3 S133RXGM 3snf uog31gtqodw3ug 31213ng aq3 20 pu1gt 390 Imammq 51gtmpgm 83 gm 133ng 39ultgt 31gtodm3ug 31213ng 1gt sagmn q1gtoxdd1gt xa33aq 31gtqmamos pu1gt uommm V 39gmg39o 033 xo ggg9391 1mg sanpm Sqx1gtau 0113 aq3 30 mtma samqu uog31gtmgxoxdd1gt apma KJGA V 39gfg g 30 mum aq3 3mm am asoddng 39aAgsuadxa 31gtqmamos 83 33 as sagws 1gt Sumxmns xo 1gtx8133ug aq3 Sump saxgnbax 3g a31gtn1gtAa o 390 52 31gt 803 138313ng Hanging sgql grew 90 8mm v0 121090 9390 812280 6390 866939 390 ma 2 392 121133 303 mowq pa31gtnq1gt3 004 1103 d WEE Z 1gtx8133ug 1gtg3tmtmdxa aq3 mpgsuog 39satqm 30 1331123 1gt mag u31gtodm3ug p1gta3su mpgsum 3138311 am 2 30 mum Kx1gtx3 qx1gt u1gt 31gt f Hanging 1gt a31gtn1gtAa 03 33mg xo aAgsuadxa 83 3g uaqm vooz gt JaqmaAoN uonmodJequI 311MB N 0313 sags qd Z smgod am up q nmq sasszzd muxom od mg anl 42 5x 0d Supmq qsgmm Mgj Mama smmmn od am up mm 368 am 11 Aquotpr m 5x 2 338 am 31 muff uogmtmsum aql lt5 c w 4 ltlt my z wdoxd maads mp 538928 qu A 117 3 0 H 7 at 58 2 smmmn od Sunmodwmg a tmx pq mx ap mq gq mag pammwum aq um muxom od p 3th 39u am ap 30 muxom od p ug Sumnsax smgod 71 up apnpug pmm 3M 13th dms Sqm 39Rxoqzzmd p Kg 03 suxgod mug saw 31 xxxmmyxm mm muxom od Supmodwmg UP pmaaps am 3g paAmdmg aq mth mm mm Spmap as 39adqs u a tmqa Jmmmxp p s mam 3219 338 am 90 pm 30 Imammq pm 30 pm 390 uaamaq Homing 319 u a tmqa am augumxa am 3g mqu mo 03 txguxmag uogqelodxatml a ue eq pazgexauaf 391 quotV WW7 8 m3 1 2137 0 0mg 0 mm7 I Mm7 Ahwdmd am GARq qu 39smmom od Sunmodwmg a tmx pq 30 saxdumxa aw asaql 058 12 058 2 117 0 0 lx 4 maqm mg71h rm7 Nod mm aq ug U633ng aq 08p um Ixonmodxmug xmug aql sgql 39awpgm am puma mmq mq axqm am 30 SWAJG XE 3SR pm 381 019 u 3808110 aq 03 8pm xmw mp q gq sg mmom od am 30 00180 am pm pamds KHRHbG aw smgod am 31 39pamadxa aq mth umq Imammq spamug am u sasgx pm smgod uogmodmm am 31gt saqsgmm putmq mm am 31219 aapoN 0 0x 0H gt W J 00 2mm am no putmq p sn saAg mm ngm f 30 aAnmgmp 381 u aq 30 aptmu mn up no a punoq mddn UP pug um am samgmmos mg 392 no spuadap 3 moq moug mp am 80823 up 1201 1 39f 30 aAnmgmp 3S u aq summ Hmf Ixonmou aql 39W Ox 120103ng mp ug sag mq 2 no spuadap 3 amqm 0 xult0d u 31uampr 3031238 31219 Immatp p q nmq 3RIAAGIHOS pagntmnb aq um Katmdamsgp am f uogmun qmoms A qmgns Km 101 moqp Sq Ami sg mmom od 391mm aql 601791391 2mg sg 8mm mnpap xgs 03 0mm 330mm aq squm 61I39I 91390 1 368 am 91390 w n m npgtAa pm GAOQR 0mm mp ug suxgod my up Sunmodwmg muxom od 391mm aq unusum am 31 3312mm aq o3 punoq sg uommmsmdm mmouxxod p as Rgmou od p mu sg mmodwmg o3 Sundmmw aw am uogmun 1218mm mum odxa aql nxmm aq 30 031211030 w mmmm Sugtmam mg mme pmgmnmx 0N uogquOdJamI ug 510113 3391 a0 Vm 6390 9 0 am 0 8390 V0 90 2 W0 am V0 I wng uogmtmsum a tmx pq aql quotamaxpmxb p aq 33mm 3g smgod mug WM um um u VO Q39O Z39O x 31gt sanpm aq Supmodwmg muxom od am mun max mxsm sgq 30 HOQRJHddR UP 8V 39mxsm anbgtm mg mng uogmtmsum mo smgod u up Sunmodwmg mmou od u 00180 mm 5pm aq um mam 03ng 39paxgnbax 8170104 W0 me0 Vm 8188390 9 V0 W210 I J 16090 AAAA mm am 30 3mm p sammua pm 2 mmd Ixogwxodwmg pmsanbax mp mm sums Immm w sgql 39mtmxo m amAaN up 3st 5mm am mmom od am 30 suxapggam am moug 03 mm mp am pm anpm pamxodmmg am 338 0 Swan sg p08 mo 31 39snogpa 3RIAAGIHOS sg mmou od Supmodwmg am Sunammmg mqwo le BIIIABN z 900 04 X37Xf79810d OOO O 900390 39smgod pamds tha Sunmodwmg Imqm smmmn od am ap q gq pgtmz pmoqs am as xxm g 319 u 338 um am uomutg Mg 319 m 8123 am 8 LO ZWF39O 9390 0189M 682260391 SWOI39O V0 6668 968160 98111 OZSZI39I AZOS39O 8390 221791139 981801 90181 812280 6390 08222281 862291 390 3991390 52 31gt pawxodmmg xm 33ng Rntmtxodxa am 30 axdumxa am 03 mumps amAaN mp Kidd am awH 39pmoqs 3g g 2 1 g 103 Hi Lxi31d 3mg gymx pmoqs umpms aql 2 1 mng amx am axdumxa 103 39333 mp 03 Immxm auo u sanpm 30 mm mag Immxm Sup u sanpm am mmmxm 03 past 8 Mimic yr 31x I xr 1 1 141quotquot g1x x x1quot 1 lt14 d x x my mnuapg aql 39smgod am up Sunmodwmg mm samm pm Sum mm 5pm Spq ummm 38m aql 39GAOQR mid pawnsnm 3M 39smgod snon num mug Sunmodwmg mag samm ummm mm aql 39GAOQR pmannsum am mmom od 319 s 11H 3219 aagoN 39smgod 3mmpr 30 mud Sunmodwmg mag gamma ummm p1th aql 39Wi samm 12ngng aq Suppm3 smnpmdax 3g 3mg Sugmmm qam mmd mm 5pm Supmodwmg mm samm ummm pumas mp 528 03 8 um i 392 mmd mp 31gt paw npgtAa 52 39 39 39 15 52 smgod am Sunmodmmg muxom od am sg amqm x W WW2 Yd 11 MEZIJ 2261 amed M22161 261 Km amed W Od rm 2 W061 833333333333 30 03338120333 33008 12 333333 803123333380 30330 33305 03312333 2333038 0333 30 31230333 0333 39033333033 33030033330 031233 33383333 033 833033123033303333 312333303353033 00380p 338333 30 8p312z1233 0333 333312 3338330 0333 312 33033033333 83333 333 53331233383338 333333333312803 0333 33303312 03 033 3123333 8333330333 33033 3930330 3303123333380 3330 003333 53312033 3000 5331203 83 03331233 38123 0333 333 30330 0333 08 60339339 83 3333833 3338 03 03331233 312333012 0333 33331203 033 0312 3330 3123 33033 003 33303312 30 339133 63133 83 33 08123 3330 333 3938033 333033 0333 333312 03331233 38033 0333 3300333033 083312330 0333 30 0p333333812333 0333 533 330338 83 30330 0333 30 03123333380 0333330 3 303123033303333 53333 33 330333 033 33 03331233 3012330 0333 312 8333gt300d 3330333333 3333803 38033 3330 333 30330 33033123033303333 0333 30 03123333380 3312 330033 033 33123333 390333123 12 333 83333033 0333 3312 833338331233330 0303033 03331283Ap12 83333333038 0312333 512333 33033333 33330d 33033123033303333 0333 30 03338 3033330 330 80333123 33031233333123 30 30333333333 0333128 0333 533383303 333333 833033123033303333 0p 03 38033 83 33 0333123 12 30 808330 0333 312033 83333033 330012338 5333030 330 8312333303353033 003803 338333 333333 8303812833 3330312 0 33330 3033333 0333 03 33033123033 03333033 0333 3330 K3312 03 Aquot3128803033 831233312 3033 83 33 83333033 0333 3312 83333123033303333 312333303353033 033312331 0333 333033 833333300 639 033 333330338 3333803 38033 3330 338333 003 33383 003803 0333 83303 381203 312 312333303353033 0333 30 003803 0333 0812030333 033 0A03du33 03 03123333380 0333 30033330 033 033335 39639 pg39 3133393 23333333300 3333333 0333 333333 8333331238 3333333300 333120 333 30333333333 38333 0333 333033 331203 033 33120 80333123 83333123033303333 0333 0A0q12 033333312330 0333 301 398033312A 83333123033303333 0333 K3330 0833 333330338 033 53312030 08 2333833 033 33123 8330331230331233x0 3123333033533 3303383312 38033 0333 80338 03331233 3303333333 31233303333 g 0333 3030 333033 8323 p3312 30063 S3123UIOUK30d 0333 033338 833033123033303333 033 33333033 3333333300 33333303 0333 333 80333123 0333 38333 0333 930 23 31gt 330331230d303333 303 0333123 3123333338 12 33030333383300 331233 033 33 033333312330 301 393333333300 12 333 80333123 83333123033303333 0333 033 33120 030333 398330331230331233330 0312 8303330 033 3933033123033303333 3312 5331203 83 3333333300 312333 333 03331233 38333 0333 K3330 gm 2 301 39830333383033 0333 833 333312 3303 312333 330 3333033 312333300 0333 30 33033123033303333 333123p1233b 12 30 3333803 0333 80338 3333333300 33333303 0333 83123303333 303330 0333 330 833033123033303333 333033 933 03 3331233 8330331230331233x0 31203333 53312333012 0312 3333333300 3333333 0333 333 8303330 033 3933033 1230d303333 3312 5331203 83 gggg393 033312A 38333 0333 K3330 312333 0033033 930 x 312 33033033333 83333333803 0333 833331233312A0 333312 83333033 833330333383033 0333 0333 3300333033 330331230d303333 31203333 12 833303 30 3333803 0333 80338 3333333300 3333333 0333 803331233 Vi 0333 30 53303 12 83 3333333300 33330008 0333 p08333303d 8V 3980333123 0333 0333 03123333303 39830330 3330 33333303 03333303 03 8330331233303120 03123330333303333 0333 333 33833 3123333003 1233330 3312 K3312 033A m n 39uum1m 3S1 1 up may papaau 1113s sanpm ammxmo mu op umnpm 113 113 m sanpm 11mg 113 as spmmpzzq g on i001 muug 113 um 03 8 mm 111 nixommn mndmm upms xgmzm 113 30 umth mm 5pm sdaag 31 as pau gsap aq um Immm m m1 mm 1315 39aspwxa 1w apm mndmm m1 pup apmopmsd 113 30 uogmnnsum m1 aAm1 3M 39Wi 076 sanpm 112131ng m1 mm pm 39 39 quota 1 1 pm 103 1 u39 39 39 11 g m 4 L1 V it 1414th x MHVOU W x m 5 mng sf m1 0 1 b3 Sunwwog 39Lquot 1lt L lt Ld 576 mm as 1W6 m5 0 26 2x m5 1 26 t tf 2 26 0 26 2x 2 26 1 26 2 26 0 16 1x 69 0 06 01 312mm up u umn1m 1131 pmz mm 13 m1 p10q 03 576 xmmn p Suganp0x3u1 38388th 3mm Sunmodwmg mg 30 ammuns xmmn m1 39ngulxou od Sumqodwmg Amt mg 30 3mm m1 mndmm o3 moq sn 1133 1 ha smgod pawwdwmg 30 3811 m1 0 mgod mom mm ppp am 8V mmwo lv 8 lPhysics 6720 Fourier Transforms November 10 2008 Fourier transforms provide information about the frequencies contained in a signal They also have many other applications in science and engineering Here we give a quick overview of the discrete Fourier transform of a real valued signal possibly the most common case We start with a signal ft For example the function could be a voltage varying with time Suppose we sample it at a series of times t for j 01 N 7 1 separated by a constant time interval At 711 7 tj For de niteness lets take tj jAt This gives us a sequence of sampled signal values f ftj A signal with a single pure frequency For a concrete example of sam pling suppose the signal is a pure cosine with angular frequency w ft Acoswt Acos27wt 1 If we sample it at times t jAt we get the sample values f Acosijt 2 The sampling domain in time is nite since we have j 01N 7 1 That is we sample over a total time interval NAt When we sample over a nite domain we need to say what we mean by a signal oscillating purely at only one frequency We require that over the length of the sampling interval we have an integer number m of complete oscillations That is wNAt N27TVAt 2mm 3 That puts a restriction on the allowable pure frequencies namely that they are also discrete w wm where Vm mNAt wm 27mm 27TmNAt 4 If we substitute this into the equation for fj we get f Acos27rjmN 5 for a signal with pure frequency Vm The At has canceled out Of course we can do a Fourier transform of any signal no matter whether its oscillations come out even Here we are just describing what we mean by a component with a single pure frequency over a nite interval in time Discrete Fourier transform The discrete Fourier transform is given by N71 H Fm Z fj6727r17mN j0 for integer values of m As we observed above7 the integer m is related to frequency through Vm mNAt 7 The angular frequency is wm 27mm We say that we have transformed a signal from the time domain to the frequency domain The same formalism applies to the Fourier transform of data f depend ing on a spatial coordinate x In this case we sample at a series of equally spaced points mj jAz7 giving us7 again the discretely sampled values fj The Fourier transform Fm is for the wavenumber km 27TmNAz7 8 or remember k 27T the wavelength Am NAzm We say that we have transformed a signal from the spatial domain to the wavenumber domain Cosine and sine transforms Another way to look at the Fourier transform is to separate it into its real and imaginary parts We use the identity 6 cosx z sin 9 So we get N71 Re Fm E f cos27rjmN j0 N71 lm Fm 7 f sin27rjmN 10 j0 lt11 Power spectrum The Fourier transform Fm is complex The power spec trum is the square of the norm of the Fourier transform Pm lleZ Re Fm2 lm Fm2 12 It gives the spectral analysis of the signal For example7 the power spectrum of a radio broadcast signal is peaked at the broadcast frequency Inverse transform The signal can be reconstructed from its Fourier trans form using 1 N71 139 fj N Z Fmem mN 13 m0 The right hand side of this equation is a linear combination of terms with coef cients Fm7 each of which oscillate at a pure frequency This can be seen from the relation EZWijmN cos27rjmN z sin27rjmN 14 To make the oscillation explicit in terms of the time variable t jAt7 notice that we can rewrite the trigonometric functions as cum coswmtj z sinumtj7 15 where we used 27ijN 27TjAth wmtj If we are reconstructing a real signal7 clearly7 only the cosine term should contribute The sine term cancels because of the symmetries discussed below Symmetries of the Fourier transform It is easy to show that the discrete Fourier transform Fm is periodic with period N That is FmN Fm 16 Also it is easy to show that for our real signal7 the complex conjugate of the Fourier transform satis es F Rm FNm 17 This means7 in particular7 that F0 is real For even N7 FNZ is also real This also means that the power spectrum has the symmetry Pm PNm 18 Counting Fourier components Now we can answer the question7 how many independent Fourier components are there We start with N real signal components in the time domain But N complex Fourier components in the frequency domain count as 2N real numbers if we consider both real and imaginary components The identity above constrains them so that there are only N independent real values involved Often programs that compute the discrete Fourier transform use the con straint to pack the Fourier transform into a vector with the same number of real components as the signal vector For example7 for N 6 the Fourier transform in the Gnu Scienti c Library packs the Fourier transform this way and for N 57 this way algorithm in Numerical Recipes by Press at al uses a slightly different order for even N It puts the two real components F0 and FN2 in the rst and second positions and then continues in the same order as above Physics 6720 Curve Fitting November 24 2008 1 Linear Least Squares This section describes the theory and practice of linear least squares also called linear regression This technique is very often used to extract the best estimates of slope and intercept parameters from a set of z y data points with error 0 Suppose we have made a series of measurements giving N data points each consisting of a triple mpg130139 where 0139 is the standard deviation of the measure ment yi Suppose further that we expect that the relationship between z and y is given by the expression y mm b 1 where m is the slope and b is the intercept Suppose nally that we want to determine the values of m and b from the data We do this by adjusting b and m to minimize the difference between the measurement and the model The difference is expressed through N 397777113397 2 X20711 2 lt2 i1 Ui We will give some justi cation for this formula in the next section Notice that this procedure tends to make the differences between the y values ie the observed values and the mzib values ie the predicted values as small as possible Note also that the weighting of the points is larger if the standard deviation is smaller The minimum of X2m 1 occurs when QXZQm 0 8X28b 0 3 The derivatives can be evaluated easily if we rewrite Eq 2 as X2mb m2 Zzz 2mem b2 17 2mny 7 Qbe Zyz 4 where the sums are N N Zm2Zzai2 2121Ug etc 5 i1 i1 In this compact notation the weighted mean value of z is i E m E 1 for example Setting the partial derivatives to zero then gives mZx2bTZ Zmy mZmbZl a V 1 V This is a simple linear system of the form M a a i where the vector at has components m and 13 the vector 81s a E 96y c 9 lt E lt gt and the symmetric matrix M is 2 13 14 lt 23 1 gt 10 The matrix M is twice the Hessian matrix for the function xa1 ma2 b The Hessian matrix is the matrix of second partial derivatives aZXZaaaa The solutions are in compact matrix form m M la 11 ltgtM 1ltgt 12 where M 1 is the inverse of the matrix M namely 1444 El 25 lt13 29022176390 290 295 In terms of components the solution is 1 7 29024217290224 m i 2 14 EzzEli m I 15 EzzEli a At the minimum the value of X2 is easily found to be X3X2m7bzy2a1393 16 If we write Am min and Ab bib k it is easy to show that 2 can be written simply as or more explicitly X2mb X3Am22m22AmAmeA1921 17 2 Maximum Likelihood and Chi Square Although the least squares method gives us the best estimate of the parameters m and b it is also very important to know how well determined these best values are In other words if we repeated the experiment many times with the same conditions what range of values of these parameters would we get To answer this question we use a maximum likelihood method We start by assuming a probability distribution for the entire set of measure ments y1y2 yN We assume that N yrziiVZcr 5 2 P 18 y17y27 72m ma is the probability that one experiment results in the set of values y1y2 yN Notice that this probability is just the product of N Gaussian normal distributions for each value of yi with mean value g and error 01 In writing the joint probability in this way we are assuming that the outcome y is not correlated with the outcome yj for i 7 j A more complicated expression would be needed if correlations were present For now we take this expression as the simplest choice The idea of maximum likelihood is to replace the ideal mean values g with the theoretically expected values gimp mm I predicted by the linear func tion The probability distribution then becomes a model conditional probability Py1 yg lemb That is assuming that the slope and intercept are m and b it gives the probability for getting the result yl yg yN in a single measurement But then we use the power of Bayes s theorem to turn it around and reinterpret it as the probability Pmbly1y2 yN that given the experimental result the linear relationship is given by the parameters m and b Dropping the list of data points we write this probability as Pm b olt exp7X2m b2 19 where N lt e gt2 x2ltmbgt Z lt20 i1 Hi The proportional to77 symbol olt is there because our reinterpretation of the prob ability in terms of m and 1 requires us to reevaluate the normalization The probability Pm b is called the likelihood function for the parameter values m and b We want to nd the values m and If that are most probable ie maximize the likelihood function Clearly this condition is equivalent to requiring that we minimize X2 and leads to the result discussed in the rst section But we also have a way to estimate the reliability of our determination of the best values m and 13 From the expressions 17 we see that Pmb cx exp7X2mb2 is similar to a normal distribution in the variables m and b That is the probability is an exponential of a quadratic in these parameters Up to now we have encountered only quadratics in a single variable Now we have a quadratic quadratic form in two variables with quadratic terms involving m2 2 and bm Once we have realized this we can use standard results to estimate the error in the best t values The standard deviation in the parameter m is determined from the formula 00 72 dm db m 7 m2Pm b 21 00 which can be shown to be simply 72 M 111 22 Likewise the error in b is just a M lm 23 The inverse M 1 is called the error matrix77 for this reason Thus the diagonal ma trix elements of M 1 give the variance of the best t parameters The off diagonal elements contain information about correlations between the best t parameter val ues as we discuss in Sec 5 Written explicitly we have 03 L2 lt24 EzzEli z a 2 2 25 Eszli z 3 Goodness of Fit The value of X2 at the minimum gives a measure of how well the straight line actually ts the data To see this consider the meaning of each term in Eq The numerator is the square of the difference between the measured value given by y and the ideal straight line value given by gimp If the discrepancy between these two values is due only to statistical scatter according to the assumed Gaussian distribution then we would expect that most of the time the numerator would be roughly the same size as the denominator In other words we would expect that the average term in the expression for X2 would be just 1 and the total would be roughly N Well actually we would expect to do just a bit better than an average of 1 since the act of optimizing the value of m and I gives us an advantage Clearly if there were only two points then we could always arrange for the straight line to go through them exactly giving xg 0 In this case the advantage is everything But the more points there are the harder it will be to get a straight line unless the data really suggest a straight line and the more we expect the value of xg to come close to N The key concept here is the number of degrees of freedom77 df It is de ned as the number of independent data points minus the number of tting parameters df Ndata 7 nparam In our case it is N7 2 So suppose we evaluate xg and get an answer that is not N What does it mean If the answer is much bigger than N we might be tempted to say that the measured points differ from a straight line by much more than would be expected from the known errors 01 Therefore the data would not justify the assumption that the ideal values lie on a straight line If the value of xg is much less than N we might be tempted to say we overestimated our errors 0139 But even though the average value ought to come out around N we must remember we are dealing with statistics here so statistical uctuations could well be the cause of the discrepancy between the actual xg and N The concepts introduced in the previous paragraphs are formalized in the anal ysis of goodness of t The question we are asking can be phrased in probabilistic terms Given a set of variables that are distributed according to the multivariate Gaussian distribution of Eq 18 what is the probability PX2dX2 that X2 com puted according to Eq 2 has a value in the range X2 X2 dxz The answer is just the integral PNX2 dy1dy2dyNPy17y27397yN 6 2 7 391 7 g1ezp2 i y 7 g2ezp2 7 i 9N 7 gNemp2 a a UJZV X A change of variables to qi 7 giempUi gives PNX2 2W7N2dq1dq2 qu expliqf qg2 q12v2l 6lgtlt2 qfiqgi q12vl 28 We can think of the integration variables as de ning components of a vector q q1q2 qN The delta function requires that the square of the length of the vector be just X2 In fact in view of the delta function the exponential can be replaced by exp7X22 and the remaining integration just gives the surface area of a sphere of radius W in an N dimensional space The nal result is PNX2 We zZWNNZFI 29 This is called the X2 distribution for N degrees of freedom df When we are tting N points to a line with two adjustable parameters we must substitute N7 2 degrees of freedom in place of N in this formula to correct for the bias we discussed above We now return to the question we asked at the beginning of this section Suppose we minimized X2 and found a value xg ls it a good t Stated in more precise statistical language we ask what is the probability that we could have gotten a value at least as large as xg as a result of chance based on the probability Eq 29 If this probability is too small ie such a large value is unlikely we might suspect that the t is bad This probability is related to an integral 0 CL PNX2 gt X3 2 PNX2dX2 30 X0 The integral gives the chance of exceeding xg with N df This proability is sometimes called the con dence level It is plotted in the graph in Fig 1 The graph is based on two numbers namely xg and the number of degrees of freedom 71 D To read the graph select the curve that corresponds to 7113 Then locate the value of xg on the top or bottom and nd where the curve crosses the vertical line corresponding to xg Read the con dence level from the axis on the left The con dence level is the probability that the observed value of xg could be equaled or exceeded by merely random uctuations If this probability is low then we could reject the straight line theory with con dence For example suppose we had 10 points 71D 8 and got xg 20 much bigger than we would have expected The graph gives a con dence level of 001 in this case That means that such a large value of X2 would be expected to occur as a result of random uctuations only 1 of the time On the other hand if we had xg 10 the con dence level would be 025 so such a large uctuation would be expected about 14 of the time The con dence level graph is based on the assumption that the probability dis tribution is Gaussian as stated If the probability distribution is different or there are correlations among the measurements we can t use this graph 4 Nonlinear Chi Square Fits The terms linear and nonlinear refer to the dependence of the tting function on the parameters and not the dependent variable x Thus if we are tting to m W W 10 4 l l l l llll l l l lll 1 2 5 10 20 50 100 Z Confidence level for fits llllllll Figure 1 Con dence level vs X2 for various 711 y a bexp7cm7 and we do not adjust c but vary only a and b7 then it is a linear t7 because of the linear dependence in a and I But if we are also varying c it becomes a nonlinear t To be a linear t7 the second and higher derivatives of the tting function with respect to all tting parameters must all be zero If the parameter dependence is nonlinear we must use a more general approach Suppose our tting function depends nonlinearly on the parameters a a1 a2 7am Then we must minimize N X2a17 03927 39 39 39 7am 7 giempa17 0127 39 39 397aNl2Uz39239 i1 41 Steepest Descent One way to minimize a general function of several variables is to use the method of steepest descent This method is based on the fact that the gradient of a function points opposite the fall line the path of steepest descent The idea is to start somewhere say at the starting vector a and then step off in the direction of steepest descent by an amount 6a iconstVX2a 32 so that the new trial parameter vector is a 6a The problem with this method is that we don t know how far to go along the line of steepest descent ie we don t know what to take for the constant Actually the problem isn t just one of choosing one constant we might want to step with a different constant for each component i We could make up an algorithm that chooses a constant or set of constants and then if after taking a step we nd that X2 increases backs up and takes a smaller step However we can do better Let s rst consider a different approach and then we ll return to this question 42 Newt0n7s Method The minimization condition can be converted into the problem of solving a nonlinear system by requiring that the derivative with respect to each of the parameters must vanish We have 8X2 N 7 em l 0 7 2 gtP 39 33 801739 02 301739 We can try to solve this nonlinear system by using Newton s method This method starts with a trial value for the parameter vector a and then seeks the vector change 6a that improves the trial value Thus we seek the solution to 8X2a 5a 8X2a m 82X2a 0 7 801739 301739 aajaak 60 34 Thus the change is found by solving the linear system Cj ZMN ak 35 k where 18X2 N 7 ghempl aghemp c i 5 i 78 36 and 182X2a N 1 8mm 81 N yiy39 lazy M39 p 7Bmp 7 1 15147 7Bmp 39 7 k 2 aajaak 2 02 aak 801739 Z 02 aajaak l Actually this linear system is the same as what we had to solve for the linear least squares problem The connection can be made more explicit by realizing that a Taylor s expansion of xa in small shifts about the vector a is just X2a6a X2a72c6a6aM6a 38 to second order in 6a Thus M is twice the Hessian matrix just as before However in the case of a linear least squares problem the Taylor series was exact and solving the linear system once led directly to the optimum value of the parameter vector In the present case solving the linear system is just a step in the Newton iteration that is supposed to lead us to the solution after a number of steps To construct the components of the vector c and the matrix M we must be able to evaluate the rst and second partial derivatives of the tting function with respect to the tting parameters Then to proceed with the Newton iteration we must solve the linear system each time we take a step The error in the ith tted parameter is found from the diagonal element of the matrix MM at the minimum of X2 a M4 39 43 Levenberg Marquardt Method The two methods we described above have problems 1 The steepest descent method has no good way to determine the length of the step 2 Newton s method is based on solving a linear system in Eq 35 The matrix to be inverted can be singular 3 Moreover unless it is started close to the minimum Newton s method sometimes leads to divergent oscillations that move away from the answer That is it overshoots and then overcompensates etc We discuss the cure to the second problem and then show how the rst and third can be used to cure each other The standard cure to the singular matrix problem is to drop the second partials in the expression for M to give N 7 7 1 exp exp M E 77 7 40 7 k 1 02 8 801739 This matrix is positive de nite so nonsingular as long as there are no more tting parameters than data points and the tting function has independent variation in the parameters eg we wouldn t be able to nd the parameters a1 and a2 separately if the tting function depends on them only as a1 a2 A brief justi cation for this modi cation is discussed by Press et al in Numerical Recipes The main argument is that the second partial derivative term contains a factor representing the departure of the measured value from the theoretical value If the parameters correspond to a reasonable t these terms should be randomly positive and negative tending to cancel so the sum should be small if the data is close to the theoretical curve Outliers can make the omitted term large but their effect on the search for the minimum should probably be discarded The problem with overshooting can be solved by a method of Levenberg and Marquardt that combines steepest descent with Newton The idea is to modify the matrix M once more to have the form Milk H Maj Maj AMm i 39 In other words we increase the diagonal matrix elements in M by a factor 1 A Now suppose A is very large Then the diagonal term dominates and the solution to Eq 35 is just 1 8X2 21 magW which is the steepest descent method but with an explicit choice for the step length in each direction On the other hand if A is very small then the method becomes more Newton like So how do we pick A We could start the Levenberg Marquardt scheme with a small value of A Then we are doing Newton If we overshoot the minimum we would get an increase in X2 Then we back up and increase A by say a factor of 10 and try again Eventually A may be very large Then the method is more like steepest descent with a step size that gets smaller each time we increase A So by continuing to increase A we are guaranteed to decrease X2 eventually If X2 decreases on the other hand we decrease A by a factor of 10 In this way the method tunes A to the circumstances As we approach the minimum we expect Newton to do better so we should nd nally that a successful search concludes with a small value of A 604739 7 How much should we fuss over a determination of the minimum of X2 As far as the determination of the parameters is concerned we should really be content with getting X2 to within less than the true minimum plus one The reason is that an increase in X2 by one corresponds roughly to a change in one standard deviation in the tting parameters so is in the noise However the estimate of the error in the parameters through the formula 39 is often quite sensitive to how close we are to the true minimum so it is a good idea to adjust the stopping criterion for the minimization in accordance with the resulting effect on the error estimate 5 Error Matrix and Correlations Let s take a closer look at the meaning of the Taylor series expansion of X2 about the minimum Suppose there are only two parameters Call them a1 and a2 Let the minimum best t values be at a and a The Taylor series expansion is then 2 1 82x2 1 82x2 82X X2 x2ltaagt ltaraigt2 8a lta27a3gt2 gag araia2a m lt43 The second partial derivatives give the matrix elements of the Hessian matrix H This matrix is twice the curvature matrix M We used the inverse of the curvature matrix previously to get the errors in the parameters 82x2 H2M 44 7 7 80180 R M 1 45 From this expression we can see that the contour lines of constant X2 are con centric ellipses centered at the best t value A one standard deviation change in the parameter values raises X2 one unit above the minimum This 10 contour is called the error ellipse If the mixed partial derivative vanishes then the major and minor axes of the ellipses are parallel to the parameter axes If not then the ellipse is rotated with respect to the axes If the ellipses are particularly eccentric and rotated as shown in Fig 2 we see that it takes a large change in parameter values along the major axis in order to get to the la contour This implies that the parameter values are strongly correlated A measure of this correlation is obtained from the elements of the inverse matrix R which is also called the error matrix 0 R12 R11R22 46 A small or zero value implies little correlation A value close to 1 corresponds to a strong correlation An example of a situation that leads to such a correlation would be a linear t to a cluster of zy points that are all distant from the y axis See Fig 3 The y intercept is especially uncertain because it is far from the points But the value of the intercept is strongly correlated with the slope since any reasonably good t must be a line that passes through the points A small change in the slope of this line requires a corresponding change in the intercept If we look at the Hessian matrix for a linear X2 t we nd that in precisely this situation the off diagonal terms are large compared with the diagonal terms so p is large H HHHHHHHHH 4i i 2 i 50 i 72 i 74 i aill lHHlHHlHHl r 4 2 O 2 4 82 Figure 2 One sigma error ellipse7 showing a strong correlation between the pa rameters a1 and a2 The ellipse is obtained from Eq 43 by setting AXZ X2 7 X2a f a3 17 ie one unit above the best value of X2 0 2000 4000 6000 8000 10000 Figure 3 A linear t producing strong correlations between slope and intercept The two lines represent sn1all departures from the best t An increase in slope can be compensated by a decrease in intercept Physics 6720 IO Methods October 30 2002 We have been using Gin and cout to handle input from the keyboard and output to the screen In these notes we discuss further useful capabilities of these standard C features how to 1 redirect output selectively 2 detect an end of le condition and check for input data conversion errors7 3 control the format number of digits7 etc of output numeric values7 and 4 read and write to a named le Now that we have introduced classes and discussed operator overloading7 it is time to reveal the true identity of Gin and cout They are actually standard instances ofthe classes istrearn and ostrearn The classes7 their instances7 and the overloading of the operators gtgt and ltlt are de ned in the ltiostreamhgt header So in these notes we will be discussing some of the useful methods of these and related classes 1 C and Unix IO Streams Unix recognizes three standard input and output streams7 called standard in put stdin7 standard output stdout 7 and standard error stderr Normally standard input is taken from the keyboard and standard output and error go to the screen With iostreamh we access these streams through the standard names CJLn7 cout7 and cerr7 respectively There is also a Clog associated with the standard error stream Unix iostream default stdin istream cin keyboard stdout ostream cout screen stderr ostream cerr screen stderr ostream clog screen Why are there so many ways of writing to the screen Suppose your program myprog generates output that a user would like to redirect to a le called outfile The redirection works like this myprog gt outfile Now suppose you wanted to your program to prompt the user for some infor mation7 or you wanted the user to see an error message Such output really should go to the screen instead of to the le outfile where the user wouldn7t notice it right away The Unix solution is to provide two output streams7 stdout and stderr7 and declare that stdout gets redirected to the le in the command above7 but stderr gets sent to the screen For example7 in your C code you would write cout ltlt datai ltlt endl cerr ltlt quotThere is an errorquot ltlt endl and the data value would then be redirected to outfile while the error message would appear on the screen Without redirection everything goes to the screen One shouldnt be misled by the name cerr which implies that it is used only for error messages It is quite proper to use the cerr stream to prompt the user for input The stream Clog could also be used for that purpose 2 Some istream methods eof and fail Here we discuss two useful istream methods 7 one for checking for an end of le and another for checking for data format errors A keyboard end of le is signaled by typing Ctrl d For a regular le it occurs when attempting to read past the last byte Here is an example of its use float aMAX while 11 lt MAX Cin gtgt an if Cineof break n Notice that we use the dot operator to access the class member Be sure to include the because eof is a member function or method It returns true if the last Cin operation hit an end of le before nding any data In this example the value of 11 after the end of le is the number of data items read into the array a A formatting error occurs when attempting to read a number containing improper characters eg attempting to create an integer from the string 314 With Cin formatting errors are silent and potentially deadly That is there is no error message and the no input value is stored So it is a good idea to check for them For this purpose the class method fail returns true when the last read operation encountered an error For example we could improve the above code by replacing the line Cin gtgt an with Cin gtgt an ifcinfail cerr ltlt quotError reading aquot ltlt 11 ltlt quotnquot 3 Formatting output Tables of output numbers are easier to read if the columns line up Here we learn about controlling numeric output format with cout It is possible7 but clumsy The ANSI C function printf is better and can be used in C We show how to use it as well 31 Formatting under ostream To some extent you can control the alignment of numbers with the special characters n for end of line and t for tab But better control of formatting in cout is accomplished through its methods and the constants of a special ios class7 which you get when you include iostream Here is an example of its use in printing a table of my values in two columns Each of the numbers is right justi ed in a eld of width seven spaces and is displayed with a decimal point and four digits past the decimal inc lude ltiostreamhgt double X10 ylO coutprecision4 Set 4 digits past the decimal coutf1agsios rightios fixed Fixed point right justified Write column header cout gtgt quotn X ynquot cout gtgt quot nquot for i O i lt 10 i coutwidth7 Set width for Xi cout ltlt xi coutwidth7 Set width for yi cout ltlt yi ltlt quotnquot Here we use the ostream methods precision flags and width to set the format The method precision takes an integer for its argument specifying the number of digits past the decimal point The method flags takes a combination of constants de ned in the ios class Here are some of them ios constant purpose right right justify left left justify xed xed point notation scienti c scienti c notation oat eld either xed or scienti c hex hexadecimal Fixed point notation displays without a power of ten Scienti c notation displays the power of ten after the letter e The default value is floatfield meaning whichever ts better To combine compatible choices simply add them as we have done Of course these methods can be called as many times as needed to vary the output format The methods precision and flags set values that apply to all subsequent output unless they are invoked with new values But the width method applies only to the next cout value after which it is reset to the default value That is why we had to call it before writing X i and again yil 32 Formatting with ANSI C printf Here is how to do the same thing with the ANSI C printf statement Output from printf goes to stdout and may be mixed with cout inc lude ltstdio hgt Write a column header Printfquotn x ynquot Print quot nn for i 01 lt10 i Printfquot74f74fnquotxiy1 Notice that we need the ANSI C header stdioh with printf Here the format is speci ed by the rst argument of the printf function The percent signs 00 in the format string introduce the format conversion speci cation There are two of them one for each value written They are taken in order reading from left to right Each output value should have a corresponding format speci cation The 74 speci es a eld width of 7 and 4 digits past the decimal point The f speci es 77 xed77 format ie not scienti c notation with powers of 10 Omitting the eld width is OK but then you cant line up the numbers The eld width is actually treated as a minimum request If the value requires more space than you allow printf will grab more space Of course the numbers won7t line up then but that is much preferable to a misleading truncation Other than format conversion speci cations any characters in the format string are copied into the output as given In this example the end of line n is such a character Here are some commonly used format conversion speci cations The speci cation must agree with the numeric type shown printf eld purpose oowdf Fixed format XXXXXX float double oowde Scienti c notation XXXXX e nnn float double oowdg Variable format float double 04qu lnteger XXXX int oos Character string Char The g format produces either f e or dtype output depending on the magnitude of the number and whether there are only zeros after the decimal point This is just like the default behavior of cout The string conversion can also take a width speci cation as in oows 4 Reading and Writing to a File The simplest way to read from or write to a le is to write your code with Gin and cout and use redirection Suppose you wanted to do both at the same time If your program is called myprog and you would like to take input from infile and write output to outfile you would run it like this myprog lt infile gt outfile But this design works only if you have only one standard input or output stream Suppose you write a program that needs values from a table in a le called tabledata and also requires prompting the user for some parameters Logically there are then two input streams The solution is to use the class ifstream7 which is derived from the class istream7 so has many of its methods The extra f reminds us that it deals with a le instead of stdin The class ofstream is used for output to a le Both of these classes are de ned in the standard header fstreamh Here are the steps required for handling a le for either input or output 1 Create an instance of ifstream or ofstream 2 Open the le Check for failure to open 3 Read from the le or write to it 4 Close the le Here is an example for input A similar approach works for output include ltfstreamhgt int main ifstream table 1 Create instance float a10 tableopenquottabledataquot 2 Open the file iftablebad Check open cerr ltlt quotCan t open tabledatanquot return 1 forint i O i lt 10 itable gtgt ai 3 Read data tableclose 4 Close the file The name of the class instance table can be whatever we want We would create a new instance of the class for each le we wanted to read Notice that we use it with gtgt just like Cin The class method Open takes the name of the le as its character string argurnent That argument can also be a character array containing the name of the le to be opened The class method Close closes the le While closing the le is not absolutely required7 it is a good habit7 since there is a usually generous upper limit to the number of les that can be kept open at once As can be guessed from the example7 the method bad checks whether the open operation succeeded It could fail if the le wasn7t there7 for example It is possible to combine steps 1 and 2 using the constructor ifstream table quottabledataquot The methods eof and fail work for the ifstream class

### BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.

### You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

## Why people love StudySoup

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "I made $350 in just two days after posting my first study guide."

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "It's a great way for students to improve their educational experience and it seemed like a product that everybody wants, so all the people participating are winning."

### Refund Policy

#### STUDYSOUP CANCELLATION POLICY

All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email support@studysoup.com

#### STUDYSOUP REFUND POLICY

StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here: support@studysoup.com

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to support@studysoup.com

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.