  if ~exist('PlDataSet'), PlDataSet=1, end

  DecayRate = 1 - 0.005;
  InitSigmaSq = 0.1;

  Mu = [0 0];
  MuSum = [0 0];
  SigmaSq = [0];
  SigmaSqSum = [0];
  Npoints = [1e-10];
  VarCounts = [0 0 0 0 0];

  Patterns = plgenerate(PlDataSet);
  MatchThreshold = 0.05;

  for trial = 1:size(Patterns,1)

    inp = Patterns(trial,:);

    if trial > 1
      distsq = sum((ones(size(Mu,1),1)*inp - Mu)'.^2);
      distsqvar = distsq ./ (2*SigmaSq);
      act = (1.0 ./ (2*pi*SigmaSq)) .* exp(-distsqvar);
      [winval,winind] = max(act);
    else
      winval = 0;
    end

    newsig = InitSigmaSq;
    newflag = (trial==1 | winval<MatchThreshold);
    if ~newflag
      countind = min(size(VarCounts,2),ceil(4*sqrt(distsqvar(winind))));
      if countind==3 & 3*VarCounts(winind,1)<VarCounts(winind,3)
        if rand(1,1)<0.05
           newflag = 1;
           newsig = SigmaSq(winind);
        end
      end
    end

    if newflag      % time to make a new bump
      if trial > 1
        winind = 1 + size(Mu,1);
      else
        winind = 1;
      end
      distsq(winind) = 0;
      MuSum(winind,:) = 0;
      SigmaSqSum(winind) = newsig;
      Npoints(winind) = 0;
      VarCounts(winind,:) = 0;
      countind = 1;
    end

    MuSum(winind,:) = MuSum(winind,:)*DecayRate + inp;
    SigmaSqSum(winind) = SigmaSqSum(winind)*DecayRate + distsq(winind);
    VarCounts = VarCounts*DecayRate;
    VarCounts(winind,countind) = VarCounts(winind,countind) + 1;
    Npoints(winind) = Npoints(winind)*DecayRate + 1;

      % Delete any bumps that aren't capturing new points
    losers = VarCounts(:,1) < 0.7;
    losers(winind) = 0;  % don't delete current winner no matter what
    MuSum = MuSum(~losers,:);
    SigmaSqSum = SigmaSqSum(~losers);
    VarCounts = VarCounts(~losers,:);
    Npoints = Npoints(~losers);

    SigmaSq = SigmaSqSum ./ Npoints;
    Mu = MuSum ./ (Npoints'*[1 1]);

    if rem(trial,10) == 0, plshow, pause(1), end
%    pause2(0.2)

  end
