# st: deconvolution

 From "Andreas Aschbacher" To statalist@hsphsun2.harvard.edu Subject st: deconvolution Date Tue, 30 Mar 2004 15:34:16 +0200 (MEST)

```Dear Nick,Stephen,fellows !

I think deconvolution works for a simple example of two lognormal components
using modified Simpson for integration:
I type:
clear
set obs 4
input y
1
0
0
0
then I type >do simpson<
~~~
capture program drop nlfaq
program nlfaq
if "`1'" == "?" {
global S_1 " m1 m2 sigma1 sigma2 "
global m1 = 1
global m2 = 1
global sigma1 = 1
global sigma2 = 1
exit
}
tempvar yh
*      global sigma1=1
*egen `X' = (_n<=34)*sum(_n*cos(\$B) + _n-1*tan(\$B))
#delimit;

gen `yh' =(0.4/12)* ((1/sqrt(2*_pi))*(1/\$sigma1)*(1/.1)* exp(-(log(.51 - .1)
- \$m1 )^2/2*(\$sigma1^2)) * (1/sqrt(2*_pi))*(1/\$sigma2)*(1/.1)* exp(-(log(.1)
- \$m2)^2/2*(\$sigma2^2))
+ 4*(1/sqrt(2*_pi))*(1/\$sigma1)*(1/.2)* exp(-(log(.51 - .2) -
\$m1 )^2/2*(\$sigma1^2)) * (1/sqrt(2*_pi))*(1/\$sigma2)*(1/.2)* exp(-(log(.2) -
\$m2)^2/2*(\$sigma2^2))
+ 2*(1/sqrt(2*_pi))*(1/\$sigma1)*(1/.3)* exp(-(log(.51 - .3) -
\$m1 )^2/2*(\$sigma1^2)) * (1/sqrt(2*_pi))*(1/\$sigma2)*(1/.3)* exp(-(log(.3) -
\$m2)^2/2*(\$sigma2^2))
+ 4*(1/sqrt(2*_pi))*(1/\$sigma1)*(1/.4)* exp(-(log(.51 - .4) -
\$m1 )^2/2*(\$sigma1^2)) * (1/sqrt(2*_pi))*(1/\$sigma2)*(1/.4)* exp(-(log(.4) -
\$m2)^2/2*(\$sigma2^2))
+ (1/sqrt(2*_pi))*(1/\$sigma1)*(1/.4)* exp(-(log(.51 - .5) - \$m1
)^2/2*(\$sigma1^2)) * (1/sqrt(2*_pi))*(1/\$sigma2)*(1/.4)* exp(-(log(.4) -
\$m2)^2/2*(\$sigma2^2))) in 1;

#delimit cr
* gen `yh'= exp(\$A) + `X' - 2 in 1
replace `yh'= \$sigma1 - \$sigma2 in 2
replace `yh'= \$m1 - \$m2 in 3
replace `yh'= \$m1 - \$sigma2 in 4

replace `1' = `yh'

end
nl faq y

di \$sigma1,\$sigma2
di \$m1,\$m2
di  \$m1,\$sigma2

I'd like to be successful for a great number of Simpson steps.
I could not bring egen-idea to work.
/Simpson rule: (b-a)/3n*{ya + 4y1+2y2 .. +2y(n-2)+4y(n-1) + yb} as integral
from a to b   /

thanks to all
aa

--
+++ 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/
```