Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2013

[Date Index] [Thread Index] [Author Index]

Search the Archive

Wald's SPRT algorithm

  • To: mathgroup at smc.vnet.net
  • Subject: [mg129801] Wald's SPRT algorithm
  • From: richardgreco at gmail.com
  • Date: Fri, 15 Feb 2013 01:56:48 -0500 (EST)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • Delivered-to: l-mathgroup@wolfram.com
  • Delivered-to: mathgroup-newout@smc.vnet.net
  • Delivered-to: mathgroup-newsend@smc.vnet.net

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]



  • Prev by Date: Re: Autocomplete popup with Mathematica 9.0.1 on Xubuntu
  • Next by Date: Re: Mathematica and Lisp
  • Previous by thread: Re: Finding Maximum without a plot
  • Next by thread: Executing external command with parameters