[R] Problem with sample()

Sun Apr 5 01:41:53 CEST 2009

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
>