# Re: st: Re: convolution/cumulative

 From "Andreas Aschbacher" To statalist@hsphsun2.harvard.edu Subject Re: st: Re: convolution/cumulative Date Thu, 25 Mar 2004 15:48:17 +0100 (MET)

```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
> michael.blasnik@verizon.net
>
>
> ----- Original Message -----
> From: "Andreas Aschbacher" <aa_surf@gmx.at>
> To: <statalist@hsphsun2.harvard.edu>
> 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/
```