#### One sided hypothesis

For a hypothesis of 50% chance of success, it is pretty obvious that only the one sided hypothesis of less than 50% success can be alternative. The hypothesis would be rejected if it takes enough trials to get the first success. This can be simulated pretty easy. It can also be calculated, but for that it seems better to estimate the chances of 1 to 5 trials and take the remainder rather then make a sum of a potential infinite number of fails before the first success. For comparison the binomial test is added. Again one sided, probability of upto one success.empprob <- function(p,n=1000) {

samps <- sample(c(FALSE,TRUE),n,replace=TRUE,p=c(1-p,p))

ww <- which(samps)

dd <- diff(c(1,ww))

sum(dd>=6)/length(dd)

}

pr <- seq(0,1,.005)

sa <- sapply(pr,empprob,n=1e4)

dbin <- function(size,prob,log=FALSE) {

if (log)

log(prob)+(size-1)*log(1-prob) else

prob*(1-prob)^(size-1)

}

pbin <- function(size,prob) {

1-sum(dbin(1:(size-1),prob))

}

sa2 <- sapply(pr,function(x) pbin(6,x))

sa3 <- pbinom(1,6,pr)

plot(x=pr,

y=sa,

type='l',

ylab='Alpha',

xlab='probability under H0: pr>p',

main='Five fail one succes')

lines(x=pr,

y=sa2,

type='l',

col='blue')

lines(x=pr,

y=sa3,

type='l',

col='red')

legend(x='topright',col=c('black','blue','red'),

lty=1,legend=c('simulation, trials until first succes',

'calculation, trials until first succes',

'Binomial'))

#### Bayes estimate

Obviously, a Bayes estimate is interesting. Using LaplacesDemon I can reuse the density which I developed previously. Note that in hindsight for one experiment it should be equivalent to a binomial. After all, there is a factor six (six permutations) in the density, but that cancels out.library('LaplacesDemon')

mon.names <- "LP"

parm.names <- c('x','y')

MyData <- list(mon.names=mon.names, parm.names=parm.names)

N <-1

Model <- function(parm, Data)

{

x <- parm[1]

y <- parm[2]

LL <- dbin(6,x,log=TRUE)+dbinom(1,6,y,log=TRUE)

LP <- sum(dbeta(parm,.5,.5,log=TRUE))+LL

Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LL,

yhat=x, parm=parm)

return(Modelout)

}

Initial.Values <- c(.5,.5)

Fit <- LaplacesDemon(Model, Data=MyData, Initial.Values)

Fit

Mean SD MCSE ESS LB Median

x 0.2184365 0.1409504 0.011367733 277.9335 0.02194684 0.1991762

y 0.2087487 0.1358410 0.009702183 314.1335 0.02145930 0.1855887

Deviance 8.8469564 1.7228087 0.124345091 341.6582 7.30631061 8.1527719

LP -4.4234782 0.8614044 0.062172546 341.6582 -6.82652015 -4.0763859

UB

x 0.5671549

y 0.5254163

Deviance 13.6530403

LP -3.6531553

#### SAS

Proc MCMC has a general distribution where one can use an own distribution. Hence it is easy to do the same Bayes estimate in SAS too. Hence I took this as a toy example to try this feature.data a;

size=6;

output;

run;

proc mcmc data=a;

parm prob .5;

prior prob ~uniform(0, 1);

ll=log(1-prob)+(size-1)*log(prob);

model size ~ dgeneral(ll);

run;

See geometric distribution, followed by negative binomial distribution.

ReplyDeleteThanks. did so.

DeleteThis comment has been removed by the author.

ReplyDelete