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

Subversion Repositories wimsdev