Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.


merit_2
Posts:
3
Registered:
4/28/12


Wald's SPRT algorithm
Posted:
Feb 15, 2013 1:56 AM


Has anyone worked with Sequential Probability Ratio Test (SPRT) in Mathematica.
I have this SPRT which looks at a BernoulliDistribution[] (see below) for successes out of trials.
My application is comparing means of two distributions and to determine if I should accept or reject a hypothesis.
I am having a hard time working out it out for my application, a PoissonDistribution and comparing the means. What I am looking to do is, if I keep taking data every minute from a sensor (sensor reading 120 cpm for exmaple), I want to be able to know if it is over a threshold or is it background using SPRT.
Here is the notebook for the Bernoulli process using BernoulliDistributions. I started with this article: "The Sequential Probability Ratio Test" Phil, Robert L. (1998)
Null; alpha = .05; beta = .05; p0 = .8; p1 = .98; k0 = alpha/(1  beta); k1 = (1  alpha)/beta; c0 = Log[k0]; c1 = Log[k1]; plt1 = Plot[{c0, c1}, {x, 0, 20}, PlotRange > {1.7 c0, 1.7*c1}, DisplayFunction > Identity]; f0[x_] := PDF[BernoulliDistribution[p0], x] f1[x_] := PDF[BernoulliDistribution[p1], x] LLR[x_] := Log[f0[x]/f1[x]] Truep = .95; sum = 0; stopflag = 0; decflag = .5; n = 0; printflag = 1; sumData = {}; Print["n", " ", "samplepoint", " ", "LLR[samplepoint]", " ", \ "sum", " ", "stopflag", " ", "decflag"] While[stopflag == 0, n = n + 1; samplepoint = Random[BernoulliDistribution[Truep]]; sum = sum + LLR[samplepoint]; AppendTo[sumData, {n, sum}]; If[sum < c0, stopflag = 1; decflag = 1;, If[sum > c1, stopflag = 1; decflag = 0;, Null];]; If[printflag > 0, Print[n, " ", samplepoint, " ", LLR[samplepoint], " ", sum, " ", stopflag, " ", decflag], Null]] plt2 = ListPlot[sumData, PlotStyle > {PointerSize[.015]}, DisplayFunction :> Identity]; Show[{plt1, plt2}, DisplayFunction > $DisplayFunction]



