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

st: deconvolution


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



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