On 04/04/2009 6:34 PM, mackas21 wrote:> 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?
Presumably h doesn't contain a vector of positive probabilities. You
didn't give us a *minimal* reproducible example, but it is at least
reproducible. Using options(error=recover) breaks at the error, and at
that point I see
Browse[1]> h
[1] 9.555555556 -1.228399425 -0.161902016 -0.086878383 -0.058343015
[6] -0.043116464 -0.033609194 -0.027105323 -0.022382548 -0.018806257
[11] -0.016012859 -0.013778429 -0.011957202 -0.010450045 -0.009187093
[16] -0.008117645 -0.007203985 -0.006417462 -0.005735913 -0.005141922
[21] -0.004621615 -0.004163804 -0.003759370 -0.003400807 -0.003081885
[26] -0.002797391 -0.002542934 -0.002314790 -0.002109784 -0.001925194
[31] -0.001758673 -0.001608193
so that's your problem. I'll leave it up to you to figure out why h
contains negatives.
Duncan Murdoch
> 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
>