Package gov.sandia.cognition.math.matrix

Examples of gov.sandia.cognition.math.matrix.Matrix


      final MultivariateGaussian priorPsi = transParticle.getPsiSS();
      final KalmanFilter kf = transParticle.getFilter();

      final int xDim = kf.getModel().getInputDimensionality();
      final Matrix H = MatrixFactory.getDefault().createMatrix(xDim, xDim * 2);
      H.setSubMatrix(0, 0, Ix);
      H.setSubMatrix(0, xDim,
          // x_{t-1}
          MatrixFactory.getDefault().createDiagonal(transParticle.getStateSample()));

      /*
       * Construct the measurement prior predictive likelihood
       * t_n (H*m^psi, d*C^psi/n)
       */
      final Vector mPriorPredMean = H.times(priorPsi.getMean());
      final Matrix mPriorPredCov = H.times(priorPsi.getCovariance()).times(H.transpose())
          .plus(Iy.scale(2d));

      // TODO FIXME
      final Matrix stPriorPredPrec = mPriorPredCov.inverse().scale(
          transParticle.getSigma2SS().getShape()/transParticle.getSigma2SS().getScale());


      MultivariateStudentTDistribution mPriorPredDist = new MultivariateStudentTDistribution(
          transParticle.getSigma2SS().getShape(),
View Full Code Here


      Vector sampledAugResponse = VectorFactory.getDefault().createVector(this.numCategories);
      for (int i = 0; i < this.numCategories; i++) {
  //        final Vector betaSample = particle.getLinearState().sample(this.random);
        final MultivariateGaussian predictivePrior = particle.getLinearState().clone();
        KalmanFilter kf = particle.getRegressionFilter(i);
        final Matrix G = kf.getModel().getA();
        predictivePrior.setMean(G.times(predictivePrior.getMean()));
        predictivePrior.setCovariance(
            G.times(predictivePrior.getCovariance()).times(G.transpose())
              .plus(kf.getModelCovariance()));
 
        // X * beta
        final double lambda = Math.exp(data.getObservedData().times(
            predictivePrior.getMean()).getElement(0));
        lambdaSum += lambda;
        final double dSampledAugResponse = -Math.log(
            -Math.log(U)/(1d+lambdaSum)
            - (data.getObservedValue().getElement(0) > 0d
                ? 0d : Math.log(this.random.nextDouble())/lambda));
        sampledAugResponse.setElement(i, dSampledAugResponse);
      }

      /*
       * Expand particle set over mixture components
       */
      for (int j = 0; j < 10; j++) {
        double categoriesTotalLogLikelihood = 0d;

        final UnivariateGaussian componentDist =
            this.evDistribution.getDistributions().get(j);
 
        for (int k = 0; k < this.numCategories; k++) {
          /*
           * TODO could avoid cloning if we didn't change the measurement covariance,
           * but instead used the componentDist explicitly.
           */
          final FruehwirthMultiParticle predictiveParticle = particle.clone();
          predictiveParticle.setAugResponseSample(sampledAugResponse);
 
          predictiveParticle.setEVcomponent(componentDist);
         
          /*
           * Update the observed data for the regression component.
           */
          predictiveParticle.getRegressionFilter(k).getModel().setC(data.getObservedData());
 
          final Matrix compVar = MatrixFactory.getDefault().copyArray(
              new double[][] {{componentDist.getVariance()}});
          predictiveParticle.getRegressionFilter(k).setMeasurementCovariance(compVar);
         
          final double logLikelihood = this.updater.computeLogLikelihood(predictiveParticle, data)
              + Math.log(this.evDistribution.getPriorWeights()[j])
View Full Code Here

      final MultivariateGaussian updatedBetaMean = particle.getPriorBeta().clone();

      final List<Double> lambdaSamples =
          new ExponentialDistribution(2).sample(random, updatedBetaMean.getInputDimensionality());
      final Matrix lambdaSamplesMatrix =
          MatrixFactory.getDenseDefault().createDiagonal(
              VectorFactory.getDenseDefault().copyValues(lambdaSamples));
     
      /*
       * To perform the following updates, we need smoothed joint samples of the previous and
       * current states (beta and the global mean), i.e. a draw from (x_t, x_{t-1} | y_t). FIXME XXX
       * The above isn't currently being done
       */
      final MultivariateGaussian priorBetaSmoothedDist = getSmoothedPriorDist(particle.getPriorBeta(),
        augResponseDist, observation, particle.getPriorPredictiveMean());
      final Vector priorBetaSmoothedSample = priorBetaSmoothedDist.sample(random);
      final MultivariateGaussian postBetaSmoothedDist = getSmoothedPostDist(particle.getPriorBeta(),
        augResponseDist, observation, particle.getPriorPredictiveMean());
      final Vector postBetaSmoothedSample = postBetaSmoothedDist.sample(random);

      final Vector priorGlobalMeanSample = particle.getPriorBeta().getMean();
      final Vector postGlobalMeanSample = particle.getPriorBeta().sample(random);

      /*
       * Perform the actual Gaussian Bayes update. FIXME This is a very poor implementation.
       */
      mvGaussianBayesUpdate(augResponseDist, priorGlobalMeanSample,
          updatedBetaMean, observation.getObservedData());

      final Vector betaMeanError = postBetaSmoothedSample.minus(priorBetaSmoothedSample);
      final ScaledInverseGammaCovDistribution updatedBetaCov = particle.getPriorBetaCov().clone();
      updateCovariancePrior(updatedBetaCov, betaMeanError);
      final Matrix betaCovSmpl = updatedBetaCov.sample(random);
      Preconditions.checkState(betaCovSmpl.getElement(0, 0) >= 0d);
      updatedBetaMean.setCovariance(lambdaSamplesMatrix.times(betaCovSmpl
          .times(updatedBetaMean.getCovariance())));

      /*
       * Now, do the above for the the global mean term.
       */
      final MultivariateGaussian updatedGlobalMean =
          particle.getPriorBeta().times(particle.getAugmentedResponseDistribution());

      mvGaussianBayesUpdate(augResponseDist,
          observation.getObservedData().times(priorBetaSmoothedSample), updatedGlobalMean,
          MatrixFactory.getDenseDefault().createIdentity(
            augResponseDist.getInputDimensionality(), augResponseDist.getInputDimensionality()));

      final Vector globalMeanError = postGlobalMeanSample.minus(priorGlobalMeanSample);
      final ScaledInverseGammaCovDistribution updatedGlobalMeanCov =
          particle.getPriorBetaCov().clone();
      updateCovariancePrior(updatedGlobalMeanCov, globalMeanError);
      final Matrix globalMeanCovSmpl = updatedGlobalMeanCov.sample(random)
          .times(updatedGlobalMean.getCovariance());
      Preconditions.checkState(globalMeanCovSmpl.getElement(0, 0) > 0d);
      updatedGlobalMean.setCovariance(globalMeanCovSmpl);

      final PolyaGammaLogitDistribution updatedParticle =
          new PolyaGammaLogitDistribution(updatedGlobalMean, updatedGlobalMeanCov);

View Full Code Here

  private MultivariateGaussian getSmoothedPostDist(MultivariateGaussian postBeta,
                                                   MultivariateGaussian augResponseDist,
                                                   ObservedValue<Vector, Matrix> observation,
                                                   Vector obsMeanAdj) {
    final Matrix C = postBeta.getCovariance();
    final Vector m = postBeta.getMean();
   
    // System design
    final Matrix F = observation.getObservedData();
    final Matrix G = MatrixFactory.getDefault().createIdentity(m.getDimensionality(), m.getDimensionality());
    final Matrix Omega = MatrixFactory.getDefault().createIdentity(m.getDimensionality(), m.getDimensionality());
   
    // Observation suff. stats
    final Matrix Sigma = augResponseDist.getCovariance();
    final Vector y = augResponseDist.getMean().minus(obsMeanAdj);
   
    final Vector a = G.times(m);
    final Matrix R = Omega;
   
    final Matrix W = F.times(Omega).times(F.transpose()).plus(Sigma);
    final Matrix FG = F.times(G);
    final Matrix A = FG.times(R).times(FG.transpose()).plus(W);
    final Matrix Wtil =
        A.transpose().solve(FG.times(R.transpose())).transpose();

    final Vector aSmooth = a.plus(Wtil.times(y.minus(FG.times(a))));
    final Matrix RSmooth =
        R.minus(Wtil.times(A).times(Wtil.transpose()));
   
    return new MultivariateGaussian(aSmooth, RSmooth);
  }
View Full Code Here

  private MultivariateGaussian getSmoothedPriorDist(MultivariateGaussian priorBeta,
                                                    MultivariateGaussian augResponseDist,
                                                    ObservedValue<Vector, Matrix> observation, Vector obsMeanAdj) {
    // Prior suff. stats
    final Matrix C = priorBeta.getCovariance();
    final Vector m = priorBeta.getMean();
   
    // System design
    final Matrix F = observation.getObservedData();
    final Matrix G = MatrixFactory.getDefault().createIdentity(m.getDimensionality(), m.getDimensionality());
    final Matrix Omega = MatrixFactory.getDefault().createIdentity(m.getDimensionality(), m.getDimensionality());
   
    // Observation suff. stats
    final Matrix Sigma = augResponseDist.getCovariance();
    final Vector y = augResponseDist.getMean().minus(obsMeanAdj);
   
    final Matrix W = F.times(Omega).times(F.transpose()).plus(Sigma);
    final Matrix FG = F.times(G);
    final Matrix A = FG.times(C).times(FG.transpose()).plus(W);
    final Matrix Wtil =
        A.transpose().solve(FG.times(C.transpose())).transpose();

    final Vector mSmooth = m.plus(Wtil.times(y.minus(FG.times(m))));
    final Matrix CSmooth =
        C.minus(Wtil.times(A).times(Wtil.transpose()));
    return new MultivariateGaussian(mSmooth, CSmooth);
  }
View Full Code Here

   * @return
   */
  private void mvGaussianBayesUpdate(MultivariateGaussian obsDist,
      Vector obsMeanComponentAdjustment, MultivariateGaussian prior, Matrix X) {
    final Vector m1 = prior.getMean();
    final Matrix c1inv = prior.getCovarianceInverse();

    final Vector m2 = obsDist.getMean().plus(obsMeanComponentAdjustment);
    final Matrix Xt = X.transpose();
    final Matrix c2inv = Xt.times(obsDist.getCovarianceInverse()).times(X);

    final Matrix Cinv = c1inv.plus(c2inv);
    final Matrix C = Cinv.inverse();

    final Vector m = C.times(c1inv.times(m1).plus(
      Xt.times(obsDist.getCovarianceInverse()).times(m2)));

    prior.setMean(m);
    prior.setCovariance(C);
  }
View Full Code Here

       * and replicate their beta sampling.
       * That would require a change here...
       */
      final MultivariateGaussian predictivePrior = particle.getLinearState().clone();
      KalmanFilter kf = particle.getRegressionFilter(particle.getCategoryId());
      final Matrix G = kf.getModel().getA();
      predictivePrior.setMean(G.times(predictivePrior.getMean()));
      predictivePrior.setCovariance(
          G.times(predictivePrior.getCovariance()).times(G.transpose())
            .plus(kf.getModelCovariance()));
      final Matrix F = kf.getModel().getC();

      final UnivariateGaussian evComponent = particle.getEVcomponent();
      final double predPriorObsMean = F.times(predictivePrior.getMean()).getElement(0)
          + evComponent.getMean();
      final double predPriorObsCov = F.times(predictivePrior.getCovariance()).times(F.transpose())
          .plus(kf.getMeasurementCovariance()).getElement(0, 0);
      particle.setPriorPredMean(predPriorObsMean);
      particle.setPriorPredCov(predPriorObsCov);

      double logLikelihood = UnivariateGaussian.PDF.logEvaluate(
View Full Code Here

      final Vector invOmegaSamples = VectorFactory.getDenseDefault().createVector(1);
      final Vector z = VectorFactory.getDenseDefault().createVector(1);
      /*
       * Now, we compute the predictive log odds, so that we can evaluate the likelihood.
       */
      final Matrix x = observation.getObservedData();
      final Vector phi = x.times(particle.getPriorBeta().getMean());
      for (int i = 0; i < y.getDimensionality(); i++) {
        final double omega = PolyaGammaLogitDistribution.sample(
          phi.getElement(i), rng);
        invOmegaSamples.setElement(i, 1d / omega);
        z.setElement(i, (y.getElement(i) - 0.5d) / omega);
      }
      final Matrix zCov = MatrixFactory.getDenseDefault().createDiagonal(invOmegaSamples);

      final MultivariateGaussian augmentedPriorPredictiveLikelihood =
          new MultivariateGaussian(z, zCov);

      particle.setAugmentedResponseDistribution(augmentedPriorPredictiveLikelihood);
View Full Code Here

    for (Entry<LogitPGParticle, ? extends Number> particleEntry : target.asMap().entrySet()) {
      final LogitPGParticle particle = particleEntry.getKey();

      final MultivariateGaussian predictivePrior = particle.getLinearState().clone();
      KalmanFilter kf = particle.getRegressionFilter();
      final Matrix G = kf.getModel().getA();
      final Matrix F = data.getObservedData();
      predictivePrior.setMean(G.times(predictivePrior.getMean()));
      predictivePrior.setCovariance(
          G.times(predictivePrior.getCovariance()).times(G.transpose())
            .plus(kf.getModelCovariance()));
      final Vector betaMean = predictivePrior.getMean();
     
      final double k_t = data.getObservedValue().getElement(0) - 1d/2d;

      final int particleCount;
      if (particleEntry.getValue() instanceof MutableDoubleCount) {
        particleCount = ((MutableDoubleCount)particleEntry.getValue()).count;
      } else {
        particleCount = 1;
      }
      for (int p = 0; p < particleCount; p++) {

        for (int j = 0; j < K; j++) {

          /*
           * Sample and compute the latent variable...
           */
          final double omega = PolyaGammaDistribution.sample(0d, this.getRandom());
          final double z_t = k_t / omega;
         
          final LogitPGParticle predictiveParticle = particle.clone();
          predictiveParticle.setPreviousParticle(particle);
          predictiveParticle.setBetaSample(betaMean);
          predictiveParticle.setLinearState(predictivePrior);
         
          predictiveParticle.setAugResponseSample(
              VectorFactory.getDefault().copyValues(z_t));

          /*
           * Update the observed data for the regression component.
           */
          predictiveParticle.getRegressionFilter().getModel().setC(F);

          final Matrix compVar = MatrixFactory.getDefault().copyArray(
              new double[][] {{1d/omega}});
          predictiveParticle.getRegressionFilter().setMeasurementCovariance(compVar);
         
          final double compPredPriorObsMean = F.times(betaMean).getElement(0) ;
          final double compPredPriorObsCov =
View Full Code Here

    final Random random = new Random(seed);
    log.info("seed=" + seed);
    ExtSamplingUtils.log.setLevel(Level.INFO);

    final double trueSigma2 = Math.pow(0.2d, 2);
    Matrix modelCovariance = MatrixFactory.getDefault().copyArray(
        new double[][] {{trueSigma2}});
    Matrix measurementCovariance = MatrixFactory.getDefault().copyArray(
        new double[][] {{trueSigma2}});

    Vector truePsi = VectorFactory.getDefault().copyValues(3d, 0.2d);

    LinearDynamicalSystem dlm = new LinearDynamicalSystem(
        MatrixFactory.getDefault().copyArray(new double[][] {{truePsi.getElement(1)}}),
        MatrixFactory.getDefault().copyArray(new double[][] {{1d}}),
        MatrixFactory.getDefault().copyArray(new double[][] {{1d}})
      );
    KalmanFilter trueKf = new KalmanFilter(dlm, modelCovariance, measurementCovariance);
    trueKf.setCurrentInput(VectorFactory.getDefault().copyValues(truePsi.getElement(0)));
   
    final double sigmaPriorMean = Math.pow(0.4, 2);
    final double sigmaPriorShape = 2d;
    final double sigmaPriorScale = sigmaPriorMean*(sigmaPriorShape - 1d);
    final InverseGammaDistribution sigmaPrior = new InverseGammaDistribution(sigmaPriorShape,
        sigmaPriorScale);
   
    final Vector phiMean = VectorFactory.getDefault().copyArray(new double[] {
        0d, 0.8d
    });
    final Matrix phiCov = MatrixFactory.getDefault().copyArray(new double[][] {
        {2d + 4d * sigmaPriorMean, 0d},
        { 0d, 4d * sigmaPriorMean}
    });
    final MultivariateGaussian phiPrior = new MultivariateGaussian(phiMean, phiCov);

View Full Code Here

TOP

Related Classes of gov.sandia.cognition.math.matrix.Matrix

Copyright © 2018 www.massapicom. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.