options crt ;
name probdist 'Computing Various Probability Distributions in TSP' ;
?
?   TSP 4.3 (September 1995)
?   Author: Bronwyn H. Hall
?
smpl  1 101 ;
trend (start=-5,step=.1) x ;
trend (start=0,step=.2) z ;
?
?   Generate standard normal and standard Laplace distributions,
?   along with normal and Laplace distributed variables. 
?
normdist 0 1 x norm ;
random nvar ;
laplace 0 1 x lap ;
set lam = 1/sqrt(2) ;
random (laplace,lambda=lam) lapvar ;
?
?   Generate T and Chi-squared distributions together with T and
?   Chi-squared distributed random variables.
?
dot(value=df) '1' '5' '10' '20' ;
  tdist df x t. ;
  random (t,df=df) t.var ;
  chisq df z chisq. ;
  genchisq df chi.var ;
enddot ;

msd nvar lapvar t1var t5var t10var t20var chi1var chi5var chi10var 
    chi20var ;
?
?  If this is run interactively, you can issue plot commands to see the
?  distributions here. 
?
? ======================================================================
proc normdist mean sdev x fx ;
local z ;
?
?    Normal distribution with mean mean and std dev. sdev.
?
  z = (x-mean)/sdev ;
  fx = (1/(sdev*sqrt(2*3.14159)))*exp(-z*z/2) ;
endproc normal ;
? ======================================================================
proc tdist df x fx ;
?
?    T distribution with degrees of freedom df.
?
  local n1 ;
  set n1 = (df+1)/2 ;
  fx = gamfn(n1)/(sqrt(df*3.14159)*gamfn(df/2)*(1+x*x/df)**n1) ;
endproc tdist ;
? ======================================================================
proc chisq df x fx ;
?
?     Chi-squared distribution with degrees of freedom df.
?
  local alpha ;
  set alpha = df/2 ;
  fx = x**(alpha-1) * exp(-x) / gamfn(alpha) ;
endproc chisq ;
? ======================================================================
proc laplace mean sdev x fx ;
?
?     Laplace distribution with mean mean and std dev. sdev.
?
local z s2 ;
  z = abs(x-mean)/sdev ;
  set s2 = sqrt(2) ; 
  fx = exp(-s2*z)/(sdev*s2) ; 
endproc laplace ;  
? ======================================================================
proc genchisq df var ;
?
?     Procedure to generate random numbers that are distributed
?  Chi-squared with df degrees of freedom.
?
  random (uniform) u ;
  cdf(inv,chisq,df=df) u var ;
endproc genchisq ;
? ===================================================================