WebSVN – wimsdev – Blame – /trunk/wims/public_html/scripts/slib/stat/poisson

Rev	Author	Line No.	Line
20	reyssat	1	!if $wims_read_parm!=slib_header
		2	!goto proc
		3	!endif
		4	slib_title=Generation of random numbers with Poisson law
		5	slib_parms=2\
		6	1,m\
		7	1,lambda
		8	slib_author=S. Lemaire and Bernadette Perrin-Riou
		9	slib_out=random data of m numbers following the Poisson law \
		10	of parameter lambda (positive number)
		11
		12
		13	slib_comment= for any nonnegative integer, P(X=i)=exp(-lambda)*(lambda^i)/i!
		14
		15	slib_example= 5,2.5
		16
		17
		18	!exit
		19
		20	:proc
		21
		22	!distribute item $wims_read_parm into slib_data, slib_lambda
		23
		24	slib_out=$empty
		25	!default slib_m=1
		26	!default slib_lambda=1
		27
		28	!!slib_ex=$[exp(-$slib_lambda)]
		29	!!slib_q=$slib_ex
		30	!!
		31	!!
		32	!!slib_U=!random 0,1 repeat $slib_m
		33	!!
		34	!!!for slib_i=1 to $slib_m
		35	!! slib_V=!item $slib_i of $slib_U
		36	!! slib_q=$slib_ex
		37	!! slib_p=$slib_ex
		38	!! slib_j=0
		39	!!
		40	!!!while $slib_V >= $slib_q
		41	!! !advance slib_j
		42	!! slib_p=$[$slib_p*$slib_a/$slib_j]
		43	!! slib_q=$[$slib_q+$slib_p]
		44	!!!endwhile
		45	!!
		46	!!slib_out=!append item $slib_j to $slib_out
		47	!!!next slib_i
		48
		49
		50	slib_data=!declosing $slib_data
		51	slib_cnt=!itemcnt slib_data
		52	!if $slib_cnt=1
		53	slib_data=1,$slib_data
		54	!endif
		55	slib_out=!exec octave poisson_rnd($slib_lambda,$slib_data)
		56
		57	slib_out=!words2items $slib_out
		58
		59
		60	slib_out =!trim $slib_out
		61

Subversion Repositories wimsdev