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/
