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

Re: st: Re: convolution/cumulative


From   "Andreas Aschbacher" <aa_surf@gmx.at>
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/



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