BeginPackage["UVW`ContSamp`"]
(* Version 2.0 08/15/97 *)

RSContinuousDistribution::usage = "
RSContinuousDistribution[functionf,a,b,n] returns a sample of size n
for the distribution with density functionf on the interval [a,b]."


RSIndependent2D::usage = "
RSIndependent2D[functionf1,a1,b1,functionf2,a2,b2,n] returns a sample of 
size n {{x1,y1},...,{xn,yn}} of a two-dimensional random vector (X,Y) , 
where X and Y are independent, X has density functionf1 on the interval 
[a1,b1], Y has density functionf2 on the interval [a2,b2]."


RSNormal2D::usage = "
RSNormal2D[Sigma1,Sigma2,rho,n] returns a sample of size n for the 
two-dimensional Gaussian vector (X,Y). The means of X and Y are null, 
their standard deviations are sigma1 and sigma2. Their correlation 
coefficient is rho."


RSUnitBall::usage = "
RSUnitBall[dim,nb] returns a sample of size n of vectors uniformly 
distributed in the unit ball of the space of dimension dim."



Begin["`Private`"]


RSContinuousDistribution[density_ , a_ , b_ , n_Integer]:=

      Block[ {nbint,delta,abscissas,ordinates,cdf,inv},

	   nbint=100;
	   delta=N[(b-a)/nbint];
	   nbint=nbint+1;
	   abscissas=N[Range[a,b,delta]];
	   ordinates=Table[ density[abscissas[[i]]] ,{i,nbint}];

	   cdf = Drop[ordinates,-1]+Drop[ordinates,1];
	   delta = delta/2;
	   cdf = cdf*delta;
	   cdf = FoldList[Plus,0.,cdf];
	   cdf = cdf/Last[cdf];
	  
	   cdf = Transpose[{cdf,abscissas}];
	   inv = Interpolation[cdf, InterpolationOrder->1];

	   Return[Table[ inv[Random[] ] , {n} ] ];
       
	   ];


RSIndependent2D[density1_,inf1_,sup1_,density2_,inf2_,sup2_,n_Integer]:=
If[N[inf1]<=N[sup1] && N[inf2]<=N[sup2],
   Return[
	  Transpose[
		    {RSContinuousDistribution[density1,N[inf1],N[sup1],n],
		     RSContinuousDistribution[density2,N[inf2],N[sup2],n]}
		   ]
	 ]
  ];




RSNormal2D[sigma1_,sigma2_,rho_,n_Integer]:=

Block[{matra={{sigma1,sigma2*rho},{0,sigma2*Sqrt[1-rho*rho]}},res={},
	      u={},s},
     Do[(
	 u = {Random[Real,{-1,1}],Random[Real,{-1,1}]};
	 s = u . u;
	 While[s>1,
	       u = {Random[Real,{-1,1}],Random[Real,{-1,1}]};
	       s = u . u;
	      ]; 
	 res = Append[res,u*Sqrt[-2Log[s]/s]]     
	 ),{n}];
       res=res.matra;
       Return[res];

     ];



RSUnitBall[dim_Integer,n_Integer]:=

Block[{res={},coord={}},
      Do[(
	  coord=Table[Random[Real,{-1,1}],{dim}]; 
	  While[coord.coord>1,
		coord=Table[Random[Real,{-1,1}],{dim}]
	       ];
	  res=Append[res,coord]
	  ),{n}];
       Return[res];

     ];

     

End[ ]


EndPackage[ ]