psim <- function(x,mu0,B=100000) {
    n <- length(x); s0 <- sd(x)
    t0 <- (mean(x)-mu0)/(s0/sqrt(n))
    vx <- matrix(rnorm(B*n,mu0,s0),B,n)
    vxbar <- rowMeans(vx)
    vs <- apply(vx,1,sd)
    vt <- (vxbar-mu0)/(vs/sqrt(n))
    mean(abs(vt)>abs(t0))
}