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. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.