Maximum Likelihood routines: 
User's manual (Version 3.0)

ADVICE: Play first with pair, simubis, simu and montecarlo_lg. See later what they are
Please send your bug report and insults to one of the author
TA, Last revised January 17, 2003

What's new in Version 3.0

It is timely to have a new release of the software.  The last one was about 6 years ago.  New routines have been introduced and a few (yes a few!) bugs found.  This user's manual does not intend to be complete, so it is very likely that there are oversight and features of the software not described.  In any case let me know of any problem you have in running the software.

What's new in Version 2.0?

The previous release (1.2) allowed you to fit Fourier spectrum (the so-called amplitude spectrum). The old release was really trimmed for the LOI use. Now together with Laurent Gizon, we made it more portable so that the routines can be used for any instrument making images in a much simpler manner. Now you only need to define the leakage matrix of the mode and the covariance matrix of the noise directly in the routines that you use. An example is given here below that you modify to your own usage. It means that the routines are now slower but are easier to implement for other usage. It might mean that you do not wish to fit moderate to high degree this way!

What's new in Version 1.2?

New features have been added to the drivers.pro. For the LOI spectra ((m,nu) diagrams) you can now also fit the amplitude spectra. WARNING!: for doing that you need to know the value and sign of the crosstalk of the mode and of the noise, in the amplitude spectra. I have attached many examples that can help you to fit (m,nu) diagram provided that YOU KNOW THE CROSSTALK. For both the single spectrum and the (m,nu) diagrams, you can choose the combination of algorithm for the convergence, with their associated tolerance.

What's new after the Davos meeting?

Now you have routines for fitting the LOI spectra, or any helioseismology instruments making images.

 In addition you can now fix parameters to given values with a simple string such as: '1000000', where you have 7 parameters and the first one is fixed. It works for any kind of spectra (Unresolved or resolved instruments).
 

How to start working with it?

First using the pointers get the following files:
  • compile.pro
  • routines.pro
  • MLEfit.pro
  • drivers.pro
  • function_rot
Then start IDL.
 
 

 You can now easily compile the procedure with compile.pro. Just use the following command:

IDL>@compile

 And that's it!

 If you need more just drop me a note

1 Main routines

1.1 MLEfit

This package is the heart of the system. These routines are for fitting any spectrum using Maximum Likelihood Estimators. The density probability in the power spectrum is assumed to a chi square with 2 degrees of freedom; in the Fourier spectra each component is assumed to be normally distributed.

 The following routines are NOT to be modified without telling the author. If you found a bug just drop me a note and I will update the package. MLEfit.pro is a package with various routines:

1.2 routines

The following routines are the ones that you should input and/or modify for your own cases. They are given for completeness as some of my demo programs need it. routines.pro contains is a package with various routines that are described in more details here. For fitting single spectrum, you should write the routine describing the spectrum such as:

Function Myspectrum,param,x
; The syntax above must NOT be different from that
;
; param array describing any parameter of the spectrum
; (linewidth,amplitude, splitting, noise, etc...)
; x is the frequency
; The following 2 statements are used to know how many
; time the function is computed
common speed,Neval
Neval=Neval+1
......your stuff......
return,yourcomputation

end

Examples are given by SPECTRA, SPECTRA02, GSPECTRA

For fitting (m,nu) diagram, you should write the routine describing the spectrum such as SPECTRALTENE (if you fit the power spectra), or CROSSTALKLG (if you fit the amplitude spectra).

1.3 drivers

There are 4 routines:

1.4 function_rot

Routine for computing rotation matrices. Needed for fitting spectra with Maximum Likelihood Estimators. It returns for any degree and angle (in rad) a 2l+1 by 2l+1 matrix using Jacobi polynomials (See Edmonds, 'Angular Momentum in Quantum Mechanics').

1.5 ct

Procedure for computing mode crosstalk for velocity measurements. Needed for fitting spectra with Maximum Likelihood Estimators. It is also needed by fitgong It returns for any degree and integration diameter (in solar raiud unit) a 2l+1 by 2l+1 matrix. It uses a routine computing Legendre polynomials.

2 Examples

Here are examples how to use drivers (FITPLOTSINGLE and FITPLOT).

2.1 multi.pro

For fitting multiplets (l=1, 2, 3, etc...). The sensitivities are given in PEAKIPHIR tune it to your own.

2.2 pair.pro

For pair l=0, l=2. The sensitivities are given in PEAKIPHIR tune it to your own. There are two routines.

2.3 montecarlo.pro

For fitting synthetic (m,nu) diagram in the Fourier spectra (which is not the power spectra!). The crosstalk of the mode is that of the LOI although you can input your own values. The covariance matrix of the noise is broken into 3 independent matrices; this is what we use for the LOI. Some algebra will be required if you want to adapt it to your own instrument. The covariance matrix will depend on the masks used for combining the pixels, on the model of the solar noise and on the type of observable (intensity or velocity). This routines works for any degree but the covariance matrices are given only for up l=3.

If you want to now more just drop me a note.

2.4 fitamp.pro

For the SOHO/LOI spectra. You have to modify the parts that load the file. It is flagged by * in the procedure. If you want to now more just drop me a note.

The crosstalk matrices are generated from our own pixel sensitivity, you might have to input them before using it for your own case. We left the LOI crosstalk matrices up to l for your own convenience. If you are clever enough you may find out that fitamp.pro looks very similar to montecarlo.pro.

2.5 fitgong.pro

Yes! You can now fit the GONG (m,nu) diagrams in the Fourier spectra. If you want to now more just drop me a note. Be warned though that the crosstalk of the noise is assumed to be the same as for the modes.

This routine is very similar to fitamp.pro. As a matter of fact the differences are rather tiny. I advise that you compute the mode crosstalk and the noise covariance yourself. This is probably the only work (big though!) that you have to do.

3 IDL operation

IDL> @compile
IDL> .run multi
% Compiled module: MULTI.
% Compiled module: SIMUBIS.
IDL> simubis,1.,0.,1.,1.,0.1,2
*******************************************

Frequency of starting points= 0.00000

*******************************************
Window size= 28.0000
plotted

Thu Jun 1 17:22:20 MET DST 1995
likelihood dlikelihood
N evaluation= 1169.00
Powell finished

Number of iterations= 10 1

Fmin= -280.78849
Thu Jun 1 17:22:39 MET DST 1995
N evaluation= 1213.00
5
For l= 2
Frequency = 0.0064327785 +/- 0.073834335

Amplitude= 1.5591076 + 0.46331763 - 0.35717614

Linewidth= 0.79765624 + 0.17840557 - 0.14579642

Splitting= 0.97748720 +/- 0.037384483

Noise= 0.10814242 +/- 0.083908733
IDL> test

give the amplitude of the l=2 mode= 1.

give the frequency of the l=2 mode= 0.

give the linewidth of the modes= 1.

give the splitting of the l=2 mode= 1.

give the background noise= 0.1

give the amplitude of the l=0 mode= 1.

give the frequency of the l=0 mode= 6.

plotted

Tue Jun 20 16:40:44 MET DST 1995
likelihood dlikelihood
N evaluation= 0.00000
Powell finished

Number of iterations= 8 1

Fmin= -174.88785
Tue Jun 20 16:41:31 MET DST 1995
N evaluation= 0.00000
7
For l=2
Frequency = 0.12832522+/- 0.083315022

Amplitude= 1.0944812+ 0.27755425 - 0.22140675

For l=0
Frequency = 5.8662282+/- 0.10459044

Amplitude= 1.1993053+ 0.40139761 - 0.30074180

Linewidth= 0.87070362+ 0.16165966 - 0.13634508

Splitting= 0.95102902+/- 0.048336407
 

Noise= 0.098365665+ 0.0091338377 - 0.0083577691