Using points responsibility information to recompute the Gaussians and their weights, maximizing
the likelihood of the mixture.
168750
protected void maximization() {
		// discard previous parameters by zeroing
		for (int i = 0; i < mixture.size; i++) {
			mixture.get(i).zero();
		}

		// compute the new mean
		for (int i = 0; i < info.size; i++) {
			PointInfo p = info.get(i);

			for (int j = 0; j < mixture.size; j++) {
				mixture.get(j).addMean(p.point,p.weights.get(j));
			}
		}
		for (int i = 0; i < mixture.size; i++) {
			GaussianGmm_F64 g = mixture.get(i);
			if( g.weight > 0 )
				CommonOps_DDRM.divide(g.mean,g.weight);
		}

		// compute new covariance
		for (int i = 0; i < info.size; i++) {
			PointInfo pp = info.get(i);
			double[] p = pp.point;

			for (int j = 0; j < mixture.size; j++) {
				GaussianGmm_F64 g = mixture.get(j);

				for (int k = 0; k < p.length; k++) {
					dx[k] = p[k]-g.mean.data[k];
				}

				mixture.get(j).addCovariance(dx, pp.weights.get(j));
			}
		}
		double totalMixtureWeight = 0;
		for (int i = 0; i < mixture.size; i++) {
			GaussianGmm_F64 g = mixture.get(i);
			if( g.weight > 0 ) {
				CommonOps_DDRM.divide(g.covariance, g.weight);
				totalMixtureWeight += g.weight;
			}
		}

		// update the weight
		for (int i = 0; i < mixture.size; i++) {
			mixture.get(i).weight /= totalMixtureWeight;
		}
	}
1165732
def _preoptimize_model(self):
        
        gaussian_model = LLT(self.data, integ=self.integ, target=self.target)
        gaussian_model.fit()
        self.latent_variables.z_list[0].start = gaussian_model.latent_variables.get_z_values()[1]
        self.latent_variables.z_list[1].start = gaussian_model.latent_variables.get_z_values()[2]

        if self.model_name2 == 't':

            def temp_function(params):
                return -np.sum(ss.t.logpdf(x=self.data, df=np.exp(params[0]), 
                    loc=np.ones(self.data.shape[0])*params[1], scale=np.exp(params[2])))

            p = optimize.minimize(temp_function,np.array([2.0, 0.0, -1.0]),method='L-BFGS-B')
            self.latent_variables.z_list[2].start = p.x[2]
            self.latent_variables.z_list[3].start = p.x[0]

        elif self.model_name2 == 'Skewt':

            def temp_function(params):
                return -np.sum(fam.Skewt.logpdf_internal(x=self.data,df=np.exp(params[0]),
                    loc=np.ones(self.data.shape[0])*params[1], scale=np.exp(params[2]),gamma=np.exp(params[3])))

            p = optimize.minimize(temp_function,np.array([2.0, 0.0, -1.0, 0.0]),method='L-BFGS-B')
            self.latent_variables.z_list[2].start = p.x[3]
            self.latent_variables.z_list[3].start = p.x[2]
            self.latent_variables.z_list[4].start = p.x[0]

        return gaussian_model.latent_variables
801886
def _M2_sparse_sym(Xvar, mask_X, Yvar, mask_Y, weights=None, column_selection=None):
    
    assert len(mask_X) == len(mask_Y), 'X and Y need to have equal sizes for symmetrization'

    if column_selection is None:
        mask_Xk = mask_X
        mask_Yk = mask_Y
        Xvark = Xvar
        Yvark = Yvar
    else:
        mask_Xk = mask_X[column_selection]
        mask_Yk = mask_Y[column_selection]
        Xvark = Xvar[:, _filter_variable_indices(mask_X, column_selection)]
        Yvark = Yvar[:, _filter_variable_indices(mask_Y, column_selection)]

    Cxxyy = np.zeros((len(mask_X), len(mask_Yk)))
    Cxxyy[np.ix_(mask_X, mask_Xk)] = _M2_dense(Xvar, Xvark, weights=weights)
    Cxxyy[np.ix_(mask_Y, mask_Yk)] += _M2_dense(Yvar, Yvark, weights=weights)

    Cxyyx = np.zeros((len(mask_X), len(mask_Yk)))
    Cxy = _M2_dense(Xvar, Yvark, weights=weights)
    Cyx = _M2_dense(Yvar, Xvark, weights=weights)
    Cxyyx[np.ix_(mask_X, mask_Yk)] = Cxy
    Cxyyx[np.ix_(mask_Y, mask_Xk)] += Cyx

    return Cxxyy, Cxyyx
573645
def plot_gaussian_projection(mu, lmbda, vecs, **kwargs):
    '''
    
    '''
    return plot_gaussian_2D(project_data(mu,vecs),project_ellipsoid(lmbda,vecs),**kwargs)
748837
def _initialize_likelihood_prior(self, positions, log_likelihoods, log_priors):
        
        func = SimpleCLFunction.from_string('''
            void compute(global mot_float_type* chain_position,
                         global mot_float_type* log_likelihood,
                         global mot_float_type* log_prior,
                         local mot_float_type* x_tmp,
                         void* data){
                
                bool is_first_work_item = get_local_id(0) == 0;
    
                if(is_first_work_item){
                    for(uint i = 0; i < ''' + str(self._nmr_params) + '''; i++){
                        x_tmp[i] = chain_position[i];
                    }
                    *log_prior = _computeLogPrior(x_tmp, data);
                }
                barrier(CLK_LOCAL_MEM_FENCE);
                    
                *log_likelihood = _computeLogLikelihood(x_tmp, data);
            }
        ''', dependencies=[self._get_log_prior_cl_func(), self._get_log_likelihood_cl_func()])

        kernel_data = {
            'chain_position': Array(positions, 'mot_float_type', mode='rw', ensure_zero_copy=True),
            'log_likelihood': Array(log_likelihoods, 'mot_float_type', mode='rw', ensure_zero_copy=True),
            'log_prior': Array(log_priors, 'mot_float_type', mode='rw', ensure_zero_copy=True),
            'x_tmp': LocalMemory('mot_float_type', self._nmr_params),
            'data': self._data
        }

        func.evaluate(kernel_data, self._nmr_problems,
                      use_local_reduction=all(env.is_gpu for env in self._cl_runtime_info.cl_environments),
                      cl_runtime_info=self._cl_runtime_info)
572414
def p(i, sample_size, weights):
        

        # Determine the initial pick values
        weight_i = weights[i]
        weights_sum = sum(weights)

        # Build a list of weights that don't contain the weight `i`. This list will
        # be used to build the possible picks before weight `i`.
        other_weights = list(weights)
        del other_weights[i]

        # Calculate the probability
        probability_of_i = 0
        for picks in range(0, sample_size):

            # Build the list of possible permutations for this pick in the sample
            permutations = list(itertools.permutations(other_weights, picks))

            # Calculate the probability for this permutation
            permutation_probabilities = []
            for permutation in permutations:

                # Calculate the probability for each pick in the permutation
                pick_probabilities = []
                pick_weight_sum = weights_sum

                for pick in permutation:
                    pick_probabilities.append(pick / pick_weight_sum)

                    # Each time we pick we update the sum of the weight the next
                    # pick is from.
                    pick_weight_sum -= pick

                # Add the probability of picking i as the last pick
                pick_probabilities += [weight_i / pick_weight_sum]

                # Multiply all the probabilities for the permutation together
                permutation_probability = reduce(
                    lambda x, y: x * y, pick_probabilities
                    )
                permutation_probabilities.append(permutation_probability)

            # Add together all the probabilities for all permutations together
            probability_of_i += sum(permutation_probabilities)

        return probability_of_i
1350223
def _symmetrize_correlograms(correlograms):
    """"""

    n_clusters, _, n_bins = correlograms.shape
    assert n_clusters == _

    # We symmetrize c[i, j, 0].
    # This is necessary because the algorithm in correlograms()
    # is sensitive to the order of identical spikes.
    correlograms[..., 0] = np.maximum(correlograms[..., 0],
                                      correlograms[..., 0].T)

    sym = correlograms[..., 1:][..., ::-1]
    sym = np.transpose(sym, (1, 0, 2))

    return np.dstack((sym, correlograms))
880909
function Gaussian(parms) {
      if (parms == null) {
        parms = {};
      }
      if (parms.pi !== void 0) {
        this.pi = parms["pi"];
        this.tau = parms["tau"];
      } else if (parms.mu !== void 0) {
        this.pi = pow(parms["sigma"], -2);
        this.tau = this.pi * parms["mu"];
      } else {
        this.pi = 0;
        this.tau = 0;
      }
      if (isNaN(this.pi) || isNaN(this.tau)) {
        throw new Error("Gaussian parms can not be NaN");
      }
    }
225582
function decision2DStumpTrain(data, labels, ix, options) {
    
    options = options || {};
    var numtries = options.numTries || 10;
    
    // choose a dimension at random and pick a best split
    var N = ix.length;
    
    var ri1= 0;
    var ri2= 1;
    if(data[0].length > 2) {
      // more than 2D data. Pick 2 random dimensions
      ri1= randi(0, data[0].length);
      ri2= randi(0, data[0].length);
      while(ri2 == ri1) ri2= randi(0, data[0].length); // must be distinct!
    }
    
    // evaluate class entropy of incoming data
    var H= entropy(labels, ix);
    var bestGain=0; 
    var bestw1, bestw2, bestthr;
    var dots= new Array(ix.length);
    for(var i=0;i<numtries;i++) {
        
        // pick random line parameters
        var alpha= randf(0, 2*Math.PI);
        var w1= Math.cos(alpha);
        var w2= Math.sin(alpha);
        
        // project data on this line and get the dot products
        for(var j=0;j<ix.length;j++) {
          dots[j]= w1*data[ix[j]][ri1] + w2*data[ix[j]][ri2];
        }
        
        // we are in a tricky situation because data dot product distribution
        // can be skewed. So we don't want to select just randomly between
        // min and max. But we also don't want to sort as that is too expensive
        // let's pick two random points and make the threshold be somewhere between them.
        // for skewed datasets, the selected points will with relatively high likelihood
        // be in the high-desnity regions, so the thresholds will make sense
        var ix1= ix[randi(0, N)];
        var ix2= ix[randi(0, N)];
        while(ix2==ix1) ix2= ix[randi(0, N)]; // enforce distinctness of ix2
        var a= Math.random();
        var dotthr= dots[ix1]*a + dots[ix2]*(1-a);
        
        // measure information gain we'd get from split with thr
        var l1=1, r1=1, lm1=1, rm1=1; //counts for Left and label 1, right and label 1, left and minus 1, right and minus 1
        for(var j=0;j<ix.length;j++) {
            if(dots[j] < dotthr) {
              if(labels[ix[j]]==1) l1++;
              else lm1++;
            } else {
              if(labels[ix[j]]==1) r1++;
              else rm1++;
            }
        }
        var t= l1+lm1; 
        l1=l1/t;
        lm1=lm1/t;
        t= r1+rm1;
        r1=r1/t;
        rm1= rm1/t;
        
        var LH= -l1*Math.log(l1) -lm1*Math.log(lm1); // left and right entropy
        var RH= -r1*Math.log(r1) -rm1*Math.log(rm1);
        
        var informationGain= H - LH - RH;
        //console.log("Considering split %f, entropy %f -> %f, %f. Gain %f", thr, H, LH, RH, informationGain);
        if(informationGain > bestGain || i === 0) {
            bestGain= informationGain;
            bestw1= w1;
            bestw2= w2;
            bestthr= dotthr;
        }
    }
    
    model= {};
    model.w1= bestw1;
    model.w2= bestw2;
    model.dotthr= bestthr;
    return model;
}
666998
def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True):
    """"""
    fig, axes = pb.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True)

    # sample inputs and outputs
    S = np.ones((20, 1))
    X = np.random.uniform(-3., 3., (20, 1))
    Y = np.sin(X) + np.random.randn(20, 1) * 0.05
    # likelihood = GPy.likelihoods.Gaussian(Y)
    Z = np.random.uniform(-3., 3., (7, 1))

    k = GPy.kern.RBF(1)
    # create simple GP Model - no input uncertainty on this one
    m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z)

    if optimize:
        m.optimize('scg', messages=1, max_iters=max_iters)

    if plot:
        m.plot(ax=axes[0])
        axes[0].set_title('no input uncertainty')
    print(m)

    # the same Model with uncertainty
    m = GPy.models.SparseGPRegression(X, Y, kernel=GPy.kern.RBF(1), Z=Z, X_variance=S)
    if optimize:
        m.optimize('scg', messages=1, max_iters=max_iters)
    if plot:
        m.plot(ax=axes[1])
        axes[1].set_title('with input uncertainty')
        fig.canvas.draw()

    print(m)
    return m
228768
def optimize_svc(data, targets):
    """"""
    def svc_crossval(expC, expGamma):
        """Wrapper of SVC cross validation.

        Notice how we transform between regular and log scale. While this
        is not technically necessary, it greatly improves the performance
        of the optimizer.
        """
        C = 10 ** expC
        gamma = 10 ** expGamma
        return svc_cv(C=C, gamma=gamma, data=data, targets=targets)

    optimizer = BayesianOptimization(
        f=svc_crossval,
        pbounds={"expC": (-3, 2), "expGamma": (-4, -1)},
        random_state=1234,
        verbose=2
    )
    optimizer.maximize(n_iter=10)

    print("Final result:", optimizer.max)
764672
def gpp(V,E):
    
    model = Model("gpp")

    x = {}
    y = {}
    for i in V:
        x[i] = model.addVar(vtype="B", name="x(%s)"%i)
    for (i,j) in E:
        y[i,j] = model.addVar(vtype="B", name="y(%s,%s)"%(i,j))

    model.addCons(quicksum(x[i] for i in V) == len(V)/2, "Partition")

    for (i,j) in E:
        model.addCons(x[i] - x[j] <= y[i,j], "Edge(%s,%s)"%(i,j))
        model.addCons(x[j] - x[i] <= y[i,j], "Edge(%s,%s)"%(j,i))

    model.setObjective(quicksum(y[i,j] for (i,j) in E), "minimize")

    model.data = x
    return model
811708
def _create_model(self, X, Y):
        

        # --- define kernel
        self.input_dim = X.shape[1]
        if self.kernel is None:
            kern = GPy.kern.RBF(self.input_dim, variance=1.)
        else:
            kern = self.kernel
            self.kernel = None

        # --- define model
        noise_var = Y.var()*0.01 if self.noise_var is None else self.noise_var
        self.model = GPy.models.GPRegression(X, Y, kernel=kern, noise_var=noise_var)

        # --- Define prior on the hyper-parameters for the kernel (for integrated acquisitions)
        self.model.kern.set_prior(GPy.priors.Gamma.from_EV(2.,4.))
        self.model.likelihood.variance.set_prior(GPy.priors.Gamma.from_EV(2.,4.))

        # --- Restrict variance if exact evaluations of the objective
        if self.exact_feval:
            self.model.Gaussian_noise.constrain_fixed(1e-6, warning=False)
        else:
            self.model.Gaussian_noise.constrain_positive(warning=False)
1367421
def solve
      set_parameters
      # Vectors holding the values used for calculating the weights:
      a = NArray.float(@n)
      b = NArray.float(@n)
      # Set an initial estimate for the probabilities of true segmentation:
      @n.times do |i|
        @weights_current[i] = @decisions[i, true].mean
      end
      # Proceed iteratively until we have converged to a local optimum:
      k = 0
      while k < max_iterations do
        # Copy weights:
        @weights_previous = @weights_current.dup
        # E-step: Estimation of the conditional expectation of the complete data log likelihood function.
        # Deriving the estimator for the unobserved true segmentation (T).
        @n.times do |i|
          voxel_decisions = @decisions[i, true]
          # Find the rater-indices for this voxel where the raters' decisions equals 1 and 0:
          positive_indices, negative_indices = (voxel_decisions.eq 1).where2
          # Determine ai:
          # Multiply by corresponding sensitivity (or 1 - sensitivity):
          a_decision1_factor = (positive_indices.length == 0 ? 1 : @p[positive_indices].prod)
          a_decision0_factor = (negative_indices.length == 0 ? 1 : (1 - @p[negative_indices]).prod)
          a[i] = @weights_previous[i] * a_decision1_factor * a_decision0_factor
          # Determine bi:
          # Multiply by corresponding specificity (or 1 - specificity):
          b_decision0_factor = (negative_indices.length == 0 ? 1 : @q[negative_indices].prod)
          b_decision1_factor = (positive_indices.length == 0 ? 1 : (1 - @q[positive_indices]).prod)
          b[i] = @weights_previous[i] * b_decision0_factor * b_decision1_factor
          # Determine Wi: (take care not to divide by zero)
          if a[i] > 0 or b[i] > 0
            @weights_current[i] = a[i] / (a[i] + b[i])
          else
            @weights_current[i] = 0
          end
        end
        # M-step: Estimation of the performance parameters by maximization.
        # Finding the values of the expert performance level parameters that maximize the conditional expectation
        # of the complete data log likelihood function (phi - p,q).
        @r.times do |j|
          voxel_decisions = @decisions[true, j]
          # Find the voxel-indices for this rater where the rater's decisions equals 1 and 0:
          positive_indices, negative_indices = (voxel_decisions.eq 1).where2
          # Determine sensitivity:
          # Sum the weights for the indices where the rater's decision equals 1:
          sum_positive = (positive_indices.length == 0 ? 0 : @weights_current[positive_indices].sum)
          @p[j] = sum_positive / @weights_current.sum
          # Determine specificity:
          # Sum the weights for the indices where the rater's decision equals 0:
          sum_negative = (negative_indices.length == 0 ? 0 : (1 - @weights_current[negative_indices]).sum)
          @q[j] = sum_negative / (1 - @weights_current).sum
        end
        # Bump our iteration index:
        k += 1
        # Abort if we have reached the local optimum: (there is no change in the sum of weights)
        if @weights_current.sum - @weights_previous.sum == 0
#puts "Iteration aborted as optimum solution was found!" if @verbose
#logger.info("Iteration aborted as optimum solution was found!")
          break
        end
      end
      # Set the true segmentation:
      @true_segmentation_vector = @weights_current.round
      # Set the weights attribute:
      @weights = @weights_current
      # As this vector doesn't make much sense to the user, it must be converted to a volume. If volume reduction has
      # previously been performed, this must be taken into account when transforming it to a volume:
      construct_segmentation_volume
      # Construct a BinVolume instance for the true segmentation and add it as a master volume to the BinMatcher instance.
      update_bin_matcher
      # Set the phi variable:
      @phi[0, true] = @p
      @phi[1, true] = @q
    end
1151312
def sparse_hessian_log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # Collect the Gaussian LDS prior terms
        J_diag, J_upper_diag = self.sparse_J_prior
        H_diag, H_upper_diag = -J_diag, -J_upper_diag

        # Collect the likelihood terms
        H_diag += self.hessian_local_log_likelihood(x)

        # Subtract a little bit to ensure negative definiteness
        H_diag -= 1e-8 * np.eye(D)

        return H_diag, H_upper_diag
975434
def _update_parameters(self):
        
        self.repr_points, self.repr_points_log = self.sampler.get_samples(self.num_repr_points, self.proposal_function, self.burn_in_steps)

        if np.any(np.isnan(self.repr_points_log)) or np.any(np.isposinf(self.repr_points_log)):
            raise RuntimeError("Sampler generated representer points with invalid log values: {}".format(self.repr_points_log))

        # Removing representer points that have 0 probability of being the minimum (corresponding to log probability being minus infinity)
        idx_to_remove = np.where(np.isneginf(self.repr_points_log))[0]
        if len(idx_to_remove) > 0:
            idx = list(set(range(self.num_repr_points)) - set(idx_to_remove))
            self.repr_points = self.repr_points[idx, :]
            self.repr_points_log = self.repr_points_log[idx]

        # We predict with the noise as we need to make sure that var is indeed positive definite.
        mu, _ = self.model.predict(self.repr_points)
        # we need a vector
        mu = np.ndarray.flatten(mu)
        var = self.model.predict_covariance(self.repr_points)
        
        self.logP, self.dlogPdMu, self.dlogPdSigma, self.dlogPdMudMu = epmgp.joint_min(mu, var, with_derivatives=True)
        # add a second dimension to the array
        self.logP = np.reshape(self.logP, (self.logP.shape[0], 1))
420538
def construct_wishart(self,v,X):
        
        self.adjust_prior(list(range(int((len(self.latent_variables.z_list)-self.ylen-(self.ylen**2-self.ylen)/2)),
            int(len(self.latent_variables.z_list)))), fam.InverseWishart(v,X))