[R] Problem with sample()

[mackas21]

Hi,
I'm having a problem using sample() within a function.

Basically I get an error reading:

Error in sample(v, 1, prob = h) : non-positive probability

Can anyone advise me as to the possible origin of this error?
Here is my code

#Discretised Gillespie algorithm function (From SMfSB, D.J. Wilkinson)

gillespied=function (N, T=100, dt=1, ...)
{
        tt=0
	n=T%/%dt
        x=N$M
        S=t(N$Post-N$Pre)
        u=nrow(S)
        v=ncol(S)
        xmat=matrix(0,ncol=u,nrow=n)
	i=1
	target=0
	repeat {
                h=N$h(x, ...)
		h0=sum(h)
		if (h0<1e-10)
			tt=1e99
		else
	                tt=tt+rexp(1,h0)
		while (tt>=target) {
			xmat[i,]=x
			i=i+1
			target=target+dt
			if (i>n)
			        return(ts(xmat,start=0,deltat=dt))
		}
                j=sample(v,1,prob=h) #################error###########
                x=x+S[,j]
        }
}

#Set initial parameters

susceptible 	= 199
infectious 	= 1
recovered 	= 0
total_pop 	= susceptible + infectious + recovered
#Time to run simulation
time 		= 1000
#Recovery rate
r 		= 1/9
#Transmission rate
B 		= 0.25 / total_pop

################Rates################

p0	= B * 0.6721747879762630000000
p1	= B * 0.0885920745659092000000
p2	= B * 0.0475394709752211000000
p3	= B * 0.0319250422239456000000
p4	= B * 0.0235931405138259000000
p5	= B * 0.0183908036724927000000
p6	= B * 0.0148319138404727000000
p7	= B * 0.0122476323793264000000
p8	= B * 0.0102907015781317000000
p9	= B * 0.0087621664064845000000
p10	= B * 0.0075394959058307200000
p11	= B * 0.0065429288357371100000
p12	= B * 0.0057182186634501300000
p13	= B * 0.0050271368777644200000
p14	= B * 0.0044419396829380100000
p15	= B * 0.0039419891495863900000
p16	= B * 0.0035116072275277300000
p17	= B * 0.0031386664156453200000
p18	= B * 0.0028136371529372800000
p19	= B * 0.0025289275577435600000
p20	= B * 0.0022784155915110200000
p21	= B * 0.0020571110286774600000
p22	= B * 0.0018609069242876300000
p23	= B * 0.0016863940046095700000
p24	= B * 0.0015307200810669000000
p25	= B * 0.0013914821958686000000
p26	= B * 0.0012666429097088400000
p27	= B * 0.0011544646324891600000
p28	= B * 0.0010534576028534100000
p29	= B * 0.0009623383079593670000
p30	= B * 0.0008799959715670280000



################Individual infection number################





v0	=0
v1	=1
v2	=2
v3	=3
v4	=4
v5	=5
v6	=6
v7	=7
v8	=8
v9	=9
v10	=10
v11	=11
v12	=12
v13	=13
v14	=14
v15	=15
v16	=16
v17	=17
v18	=18
v19	=19
v20	=20
v21	=21
v22	=22
v23	=23
v24	=24
v25	=25
v26	=26
v27	=27
v28	=28
v29	=29
v30	=30




#Set conditions of Petri Net

N=list()
N$M=c(susceptible, infectious, recovered)


################Pre-matrix################

N$Pre=matrix(c(0, 1, 0, v0, 1 , 0 , v1, 1 , 0 , v2, 1 , 0 , v3, 1 , 0 , v4,
1 , 0 , v5, 1 , 0 , v6, 1 , 0 , v7, 1 , 0 , v8, 1 , 0 , v9, 1 , 0 , v10, 1 ,
0 , v11, 1 , 0 , v12, 1 , 0 , v13, 1 , 0 , v14, 1 , 0 , v15, 1 , 0 , v16, 1
, 0 , v17, 1 , 0 , v18, 1 , 0 , v19, 1 , 0 , v20, 1 , 0 , v21, 1 , 0 , v22,
1 , 0 , v23, 1 , 0 , v24, 1 , 0 , v25, 1 , 0 , v26, 1 , 0 , v27, 1 , 0 ,
v28, 1 , 0 , v29, 1 , 0 , v30, 1 , 0), ncol = 3, byrow = TRUE)




################Post-matrix################

N$Post=matrix(c(0, 0, 1, 0, v0+1, 0 , 0, v1+1, 0 , 0, v2+1, 0 , 0, v3+1, 0 ,
0, v4+1, 0 , 0, v5+1, 0 , 0, v6+1, 0 , 0, v7+1, 0 , 0, v8+1, 0 , 0, v9+1, 0
, 0, v10+1, 0 , 0, v11+1, 0 , 0, v12+1, 0 , 0, v13+1, 0 , 0, v14+1, 0 , 0,
v15+1, 0 , 0, v16+1, 0 , 0, v17+1, 0 , 0, v18+1, 0 , 0, v19+1, 0 , 0, v20+1,
0 , 0, v21+1, 0 , 0, v22+1, 0 , 0, v23+1, 0 , 0, v24+1, 0 , 0, v25+1, 0 , 0,
v26+1, 0 , 0, v27+1, 0 , 0, v28+1, 0 , 0, v29+1, 0 , 0, v30+1, 0), ncol = 3,
byrow = TRUE)




################Rate function################




N$h=function(y, th=c(r, p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11,
p12, p13, p14, p15, p16, p17, p18, p19, p20, p21, p22, p23, p24, p25, p26,
p27, p28, p29, p30))






{
 return(c(th[1]*y[2], th[2]*y[1]*y[2], th[3]*y[1]*y[2], th[4]*y[1]*y[2],
th[5]*y[1]*y[2], th[6]*y[1]*y[2], th[7]*y[1]*y[2], th[8]*y[1]*y[2],
th[9]*y[1]*y[2], th[10]*y[1]*y[2], th[11]*y[1]*y[2], th[12]*y[1]*y[2],
th[13]*y[1]*y[2], th[14]*y[1]*y[2], th[15]*y[1]*y[2], th[16]*y[1]*y[2],
th[17]*y[1]*y[2], th[18]*y[1]*y[2], th[19]*y[1]*y[2], th[20]*y[1]*y[2],
th[21]*y[1]*y[2], th[22]*y[1]*y[2], th[23]*y[1]*y[2], th[24]*y[1]*y[2],
th[25]*y[1]*y[2], th[26]*y[1]*y[2], th[27]*y[1]*y[2], th[28]*y[1]*y[2],
th[29]*y[1]*y[2], th[30]*y[1]*y[2], th[31]*y[1]*y[2], th[32]*y[1]*y[2]))
}


#Run Gillespie function

out = gillespied(N,T=time,dt=1)

simulations = 10000

for (x in (1:simulations))
{
 out = gillespied(N,T=time,dt=1)
}


Thanks in advance

Paul

-- 
View this message in context: http://www.nabble.com/Problem-with-sample%28%29-tp22888525p22888525.html
Sent from the R help mailing list archive at Nabble.com.