makeps3prob2 = function(ntrials)
{
    #In ps3prob2 the density has range [0,5] and is a mixture of gaussians
    #f(x) = b1*dnorm(x,m1,s1) + b2*dnorm(x,m2,s2) + .3*delta(5)
    #b1 and b2 are constants giving the probability is with each gaussian
    #The cdf F(x) = b1*(dnorm(x, m1,s1) - dnorm(0,m1,x1)) + b2*(dnorm(x,m2,s2) - dnorm(0,m2.s2)) + b3*step(x,5)  for x in [0,5]  --step(x,5) = unit step at 5..
    # The code generates random values by the following algorithm
    #   Pick p in Unif(0,1)
    #   Return x = F^(-1)(p)
    # The details are x = seq(0,5,.01) \approx [0,5]
    # FX = F(x) \approx F([0,5])
    # FI \approx F^(-1)([0,1])  FI[j] = largest x with F(x) < j/resolution
    # Generate uniform random samples on 1:resolution, FI[samples] are our random values

    set.seed(7)
    resolution = 10000
    m1 = .5
    s1 = .3
    a1 = pnorm(5,m1,s1) - pnorm(0,m1,s1)
    b1 = .5/a1
    m2 = 5
    s2 = 2.5/3
    a2 = pnorm(5,m2,s2) - pnorm(0,m2,s2)
    b2 = .3/a2
    x = seq(0,5,.01)

    #This plots the continuous part of the density exactly
    y = b1*dnorm(x,m1,s1) + b2*dnorm(x,m2,s2)
    plot(x,y, type='l')

    FX = b1*(pnorm(x,m1,s1)-pnorm(0,m1,s1)) + b2*(pnorm(x,m2,s2)-pnorm(0,m2,s2))
    FX[length(FX)] = 1  #add the remaining probability to F(5)
    currInd = 1
    FI = rep(0,resolution)
    for (j in 1:resolution)
    {
        p = j/resolution
        while(TRUE)
        {
            if (FX[currInd] >= p)
            {
                FI[j] = x[currInd]
                break
            }
            else
                currInd = currInd + 1
        }
    }
    s = sample(1:resolution,ntrials,replace=TRUE)
    return(FI[s])
}

# Inside R studio we generated the data ps3prob2.r with the following code
#x = makeps3prob2(5000)
#n = length(x)
#outname = 'ps3prob2data-b.r'
#pre = paste("getprob2data = function()\n{\n x = c(",x[1])
#post = paste(x[n], ")\nreturn(x)\n}")
#cat(pre, x[2:(n-1)], post,file=outname,sep=',')