Yes Michael you're right.
I also improved the case `j' = x1 /because of log / ::
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clear
set more off
set obs 200
gen d1 =.
scalar m1 = 3.91
scalar m2 = 4.6
scalar sigma1 = 0.336
scalar sigma2 = 1.609
gen x = _n
gen x1 = x*.1
forvalues i = 1(1)200 {
local j = `i' * .1
if `i' == 1 {
local j = .11
di `j'
}
gen double y =(1/sqrt(2*3.14))*(1/sigma1)*(1/x1)* exp(-(log((`j') - x1) - m1
)^2/2*(sigma1^2)) * (1/sqrt(2*3.14)) *(1/sigma2)*(1/x1)* exp(-(log((`j') -
x1) - m2)^2/2*(sigma2^2))
replace y = 0 if x1 > `j'
replace y =(1/sqrt(2*3.14))*(1/sigma1)*(1/x1)* exp(-(log((`j'+ 0.01) - x1) -
m1 )^2/2*(sigma1^2)) * (1/sqrt(2*3.14)) *(1/sigma2)*(1/x1)* exp(-(log((`j'+
0.01) - x1) - m2)^2/2*(sigma2^2)) if `j' == x1
di `j'
set more on
integ y x1 ,gen(Sy)
replace d1 = Sy if [_n] == `i'
drop Sy y
}
integ d1 x1,gen(Sy_1)
*d1 gegen x1 liefert gesuchte gefaltete Wahrscheinlichkeitsverteilung:
twoway line d1 x1
graph save "C:\DATA\convolution.gph",replace
twoway line Sy_1 x1
graph save "C:\DATA\cumulative.gph",replace
graph combine convolution.gph cumulative.gph,title(" convolution mit Stata
") note(" Bildung der Kumulativen ")
graph save "C:\DATA\Convolution_mit_Stata.gph",replace
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Now I am looking for a way to go the way back,finding the two basic curves
.............
andreas aschbacher
Just a couple of comments --
>
> You use some strange syntax which appears to be allowed by Stata, but is
> certainly not standard :
>
> gen x = [_n]
>
> is usually written: gen x=_n
>
> and
>
> replace d1 = Sy if [_n] == `i'
>
> would be better written: replace d1 = Sy in `i'
>
> Also, why do you generate x and x1 within the loop and then drop them at
> the
> end of the loop just to recreate again at the top of the loop and finally
> after the loop? -- they are always the same, aren't they?
>
>
> Michael Blasnik
> [email protected]
>
>
> ----- Original Message -----
> From: "Andreas Aschbacher" <[email protected]>
> To: <[email protected]>
> Sent: Thursday, March 25, 2004 5:55 AM
> Subject: st: convolution/cumulative
>
>
> > Dear fellows !
> >
> > maybe the following program is useful for someone:
> > / computing convolution-integral and cumulative-function /
> > I seek for possibilities to show that dosis-distributions in population
> > are composed via convolution of two functions(lognormal ?)
> > one function represents background - that is radioactivity which is
> > always around us and not-naturally radioactivity which is for example
> > radioactivity in a hospital.
> > physically both allotments are the sum - but !!!
> probability-distribution
> in
> > human-population
> > is - as we all suppose - composed via convolution.(Bayes-statistic)
> >
> > ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~+
> > clear
> > set more off
> > set obs 100
> > gen d1 =.
> > scalar m1 = 3
> > scalar m2 = 3.5
> > forvalues i = 1(1)100 {
> > gen x = [_n]
> > gen x1 = x*.1
> >
> > local j = `i' * .1
> >
> > if `i' == 1 {
> > local j = .11
> > di `j'
> > }
> > di `j'
> > di `j'
> >
> > gen double y = exp(-(log((`j') - x1) - m1 )^2/2) *exp(-(log((`j') - x1)
> -
> > m2)^2/2)
> > /standard-deviation = 1 for both /
> > replace y = 0 if x1 > `j'
> > di `j'
> > integ y x1 ,gen(Sy)
> > replace d1 = Sy if [_n] == `i'
> > drop Sy y x x1
> > }
> > gen x = [_n]
> > gen x1 = x*.1
> >
> > integ d1 x1,gen(Sy_1)
> > twoway line d1 x1||line Sy_1 x1
> >
> > ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> > greetings from austria
> > andreas aschbacher
>
>
> *
> * 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/
>
--
+++ NEU bei GMX und erstmalig in Deutschland: T�V-gepr�fter Virenschutz +++
100% Virenerkennung nach Wildlist. Infos: http://www.gmx.net/virenschutz
*
* 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/
