Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: LevenbergMarquardt


From   "Andreas Aschbacher" <[email protected]>
To   [email protected]
Subject   st: LevenbergMarquardt
Date   Tue, 20 Jan 2004 15:39:39 +0100 (MET)

Dear fellows :
In austrian centre of physics we use Stata to do
LevenbergMarquardt(nl-command)since about
four months,about 400 measurements were evaluated.
we consider small (not infinitesimal) dosis-intervals and the number of
persons which is exposed
to these intervals - the following little program shows how to count the
number of persons in certain
intervals :

/pos is ln(dosis-value) because then it seems to be Gau�-distribution
without logarithmusnaturalis it is lognormal distribution -
it is easier to reach convergence with Gau�-formula in
most cases /

local i = 3
while `i' <= 20 {
        count if pos >= `i' & pos  < `i' + .05  
        display
        local i = `i' + .05
}
*  3           -            3.05                ->  3.025
*  3.05      -            3.1                  ->  3.075
*  3.1        -            3.15                ->  3.125
*  3.15      -            3.2                  ->  3.175
*  3.2        -            3.25                ->  3.225
*  3.3        -            3.35                ->  3.325

this program counts and you get  also zero if there is no value in interval
with these values I can obtain the following ascii-text-file :

Ascii-Text-file :  it is called M012.txt below
x              y                       x...dosis-value in radioactive
nmeasurement    y...number of counts 
3.025	1                       ->for example  :
3.075	2                       two counts in interval 3.05  to 3.1
3.125	4
3.175	7
3.225	6
3.275	5
3.325	10
3.375	13
3.425	22
3.475	13
3.525	24
3.575	41
3.625	41
3.675	60
3.725	67
3.775	81
3.825	98
3.875	143
3.925	162
3.975	181
4.025	241
4.075	268
4.125	301
...
->
now import ascii-file and start LevenbergMarquardt(nl-command) with do
zweigau�_1
->
. clear

. insheet x y using "C:\DATA\Ergebnisse_5\M012.txt", tab
(2 vars, 161 obs)

. do zweiGau�_1

. capture program drop nlexample

. program nlexample
  1. version 8
  2. if "`1'" == "?" {
  3. global S_1 "A1 M1 S1 A2 M2 S2"
  4. global A1 = 446.12                
  5. global M1 = 4.34
  6. global S1 = 0.426
  7. global A2 = 28.04
  8. global M2 = 5.2
  9. global S2 = 0.94   
10. exit
 11. }
 12. replace `1' = $A1*exp(-((x-$M1)/$S1)^2) + $A2*exp(-((x-$M2)/$S2)^2) 
 13. end
/......
*I try to find good start values with :
*gen y1 = A1*exp(-((x-M1)/S1)^2) + A2*exp(-((x-M2)/S2)^2) 
*with :integ y x and :integ y1 x    I can find good start-values very easy:
* difference of the two areas should be very small i.e. less than 0.5%
*these functions(sum of two Gau�ians) are very sensitive in convergence ->
*if I use 5.19 instead of 5.2 for M2 the algorithm does not converge !!
                                                                            
                       ...........      /
. nl example y x,level(95) leave eps(1e-o5) trace iterate(15000)
delta(4e-10)
(obs = 161)

Iteration 0:  
           A1 =    446.12
           M1 =      4.34
           S1 =      .426
           A2 =     28.04
           M2 =       5.2
           S2 =       .94
                                         residual SS =  15304.15

Iteration 1:  
           A1 =  428.9222
           M1 =  4.343089
           S1 =  .4179389
           A2 =  30.91823
           M2 =  4.843218
           S2 =  1.243443
                                         residual SS =  13197.31

................................................................

Iteration 15:  
           A1 =  392.4705
           M1 =  4.364213
           S1 =  .3863484
           A2 =  68.03155
           M2 =  4.494176
           S2 =  .9563224
                                         residual SS =  6895.428

Iteration 16:  
           A1 =  392.4681
           M1 =  4.364213
           S1 =  .3863467
           A2 =  68.03416
           M2 =  4.494168
           S2 =  .9563062
                                         residual SS =  6895.428

Iteration 17:  
           A1 =  392.4666
           M1 =  4.364213
           S1 =  .3863457
           A2 =  68.03578
           M2 =  4.494164
           S2 =  .9562961
                                         residual SS =  6895.428


      Source |       SS       df       MS                     Number of obs
=       161
-------------+------------------------------               F(  6,   155) =  
8504.37
       Model |  2269985.57     6  378330.929         Prob > F      =   
0.0000
    Residual |  6895.42818   155  44.4866334       R-squared     =    0.9970
-------------+------------------------------               Adj R-squared =  
 0.9969
       Total |     2276881   161   14142.118           Root MSE      =  
6.66983
                                                                       Res.
dev.     = 
1061.809
(example)
------------------------------------------------------------------------------
           y |      Coef.      Std. Err.          t         P>|t|      [95%
Conf.
Interval]
-------------+----------------------------------------------------------------
          A1 |   392.4659   8.621851      45.52    0.000     375.4344   
409.4974
          M1 |   4.364213   .0027606  1580.89    0.000      4.358759   
4.369666
          S1 |   .3863452   .0065244      59.22    0.000     .3734571   
.3992334
          A2 |   68.03646    8.93421        7.62     0.000     50.38793   
85.68498
          M2 |   4.494162   .0286298    156.97     0.000     4.437607   
4.550717
          S2 |   .9562918   .0576867      16.58     0.000    .8423382   
1.070245
------------------------------------------------------------------------------
 (SEs, P values, CIs, and correlations are asymptotic approximations)

. predict myvar
(option yhat assumed; fitted values)

. twoway line y x || line myvar x
. integ y x
. integ myvar x  

Now we make background correctur based on these values.
The great hump is assumed to be background.
Just imagine y(x) as a snappy curve mostly composed of two or three humps
and myvar(x)(fitted) as clogged curve.
this method is to be enhanced furthermore.

Because of my first course in Stata on 23th january(it was not possible
earlier)
I had to ask many "stupid" questions ...
I am especially indebted  to Nick Cox;after netcourse I will create a
website
                                                                      
andreas

-- 
+++ GMX - die erste Adresse f�r Mail, Message, More +++
Bis 31.1.: TopMail + Digicam f�r nur 29 EUR http://www.gmx.net/topmail

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2024 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index