Two-way mixed ANOVA.
325201
def _mixed_anova(self, dv=None, between=None, within=None, subject=None,
                 correction=False, export_filename=None):
    """"""
    aov = mixed_anova(data=self, dv=dv, between=between, within=within,
                      subject=subject, correction=correction,
                      export_filename=export_filename)
    return aov
881816
def _psi2computations(dL_dpsi2, variance, Z, mu, S):
    
    # here are the "statistics" for psi1 and psi2
    # Produced intermediate results:
    # _psi2_dvariance      Q
    # _psi2_dZ             MxQ
    # _psi2_dmu            NxQ
    # _psi2_dS             NxQ

    variance2 = np.square(variance)
    common_sum = np.dot(mu,(variance*Z).T)
    if len(dL_dpsi2.shape)==2:
        Z_expect = (np.dot(dL_dpsi2,Z)*Z).sum(axis=0)
        dL_dpsi2T = dL_dpsi2+dL_dpsi2.T
        common_expect = np.dot(common_sum,np.dot(dL_dpsi2T,Z))
        Z2_expect = np.inner(common_sum,dL_dpsi2T)
        Z1_expect = np.dot(dL_dpsi2T,Z)
    
        dL_dvar = 2.*S.sum(axis=0)*variance*Z_expect+(common_expect*mu).sum(axis=0)
    
        dL_dmu = common_expect*variance
    
        dL_dS = np.empty(S.shape)
        dL_dS[:] = Z_expect*variance2
    
        dL_dZ = variance2*S.sum(axis=0)*Z1_expect+np.dot(Z2_expect.T,variance*mu)
    else:
        N,M,Q = mu.shape[0],Z.shape[0],mu.shape[1]
        dL_dpsi2_ = dL_dpsi2.sum(axis=0)
        Z_expect = (np.dot(dL_dpsi2.reshape(N*M,M),Z).reshape(N,M,Q)*Z[None,:,:]).sum(axis=1)
        dL_dpsi2T = dL_dpsi2_+dL_dpsi2_.T
        dL_dpsi2T_ = dL_dpsi2+np.swapaxes(dL_dpsi2, 1, 2)
        common_expect = np.dot(common_sum,np.dot(dL_dpsi2T,Z))
        common_expect_ = (common_sum[:,:,None]*np.dot(dL_dpsi2T_.reshape(N*M,M),Z).reshape(N,M,Q)).sum(axis=1)
        Z2_expect = (common_sum[:,:,None]*dL_dpsi2T_).sum(axis=1)
        Z1_expect = np.dot(dL_dpsi2T_.reshape(N*M,M),Z).reshape(N,M,Q)
    
        dL_dvar = 2.*variance*(S*Z_expect).sum(axis=0)+(common_expect_*mu).sum(axis=0)
    
        dL_dmu = common_expect_*variance
    
        dL_dS = np.empty(S.shape)
        dL_dS[:] = variance2* Z_expect
    
        dL_dZ = variance2*(S[:,None,:]*Z1_expect).sum(axis=0)+np.dot(Z2_expect.T,variance*mu)

    return dL_dvar, dL_dmu, dL_dS, dL_dZ
80098
function mix (r, s, or, cl) {
  if (!s || !r) {
    return r;
  }

  var i = 0,
    c, len;

  or = or || arguments.length === 2;

  if (cl && (len = cl.length)) {
    for (; i < len; i++) {
      c = cl[i];
      if ((c in s) && (or || !(c in r))) {
        r[c] = s[c];
      }
    }
  } else {
    for (c in s) {
      if (or || !(c in r)) {
        r[c] = s[c];
      }
    }
  }
  return r;
}
389421
def evaluate(x, y, constant, amplitude, x_mean, y_mean, x_stddev,
                 y_stddev, theta):
        """"""

        model = Const2D(constant)(x, y) + Gaussian2D(amplitude, x_mean,
                                                     y_mean, x_stddev,
                                                     y_stddev, theta)(x, y)
        return model
393534
def covariance2(xs, ys)
      raise "Length xs and ys must be equal" unless xs.length == ys.length

      ev_x, ev_y = mean(xs), mean(ys)
      xs.zip(ys)
        .map{|x,y| (x-ev_x) * (y-ev_y)}
        .inject(0) {|sum, x| sum += x} / xs.length
    end
884310
def _twofilter_smoothing_ON(self, t, ti, info, phi, lwinfo, return_ess,
                               modif_forward, modif_info):
        
        if modif_info is not None:
            lwinfo += modif_info
        Winfo = rs.exp_and_normalise(lwinfo)
        I = rs.multinomial(Winfo)
        if modif_forward is not None:
            lw = self.wgt[t].lw + modif_forward
            W = rs.exp_and_normalise(lw)
        else:
            W = self.wgt[t].W
        J = rs.multinomial(W)
        log_omega = self.model.logpt(t + 1, self.X[t][J], info.hist.X[ti][I])
        if modif_forward is not None:
            log_omega -= modif_forward[J]
        if modif_info is not None:
            log_omega -= modif_info[I]
        Om = rs.exp_and_normalise(log_omega)
        est = np.average(phi(self.X[t][J], info.hist.X[ti][I]), axis=0,
                         weights=Om)
        if return_ess:
            return (est, 1. / np.sum(Om**2))
        else:
            return est
1164754
def K(self, X, X2, target):
        """"""
        if X2 is None: X2 = X
       # i1 = X[:,1]
       # i2 = X2[:,1]
       # X = X[:,0].reshape(-1,1)
       # X2 = X2[:,0].reshape(-1,1)
        dist = np.abs(X - X2.T)
        
        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU

        k1 = np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
        k2 = (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2 
        k3 = np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )

        np.add(self.varianceU*self.varianceY*(k1+k2+k3), target, target)
18825
def Kdiag(self, X, target):
        """"""
        self._K_diag_computations(X)
        target+= self.variance*self._K_diag_dvar
864897
private static DataTable2D bivariateMatrix(Dataframe dataSet, BivariateType type) {
        DataTable2D bivariateMatrix = new DataTable2D();
        
        //extract values of first variable
        Map<Object, TypeInference.DataType> columnTypes = dataSet.getXDataTypes();
        Object[] allVariables = columnTypes.keySet().toArray();
        int numberOfVariables = allVariables.length;
        
        TransposeDataList transposeDataList;
        for(int i=0;i<numberOfVariables;++i) {
            Object variable0 = allVariables[i];
            if(columnTypes.get(variable0)==TypeInference.DataType.CATEGORICAL) {
                continue;
            }
            
            transposeDataList = new TransposeDataList();
            
            //extract values of first variable
            transposeDataList.put(0, dataSet.getXColumn(variable0));
            
            for(int j=i;j<numberOfVariables;++j) {
                Object variable1 = allVariables[j];
                if(columnTypes.get(variable1)==TypeInference.DataType.CATEGORICAL) {
                    continue;
                }
            
                transposeDataList.put(1, dataSet.getXColumn(variable1));
                
                double value = 0.0;
                if(type==BivariateType.COVARIANCE) {
                    value = Descriptives.covariance(transposeDataList, true);
                }
                else if(type==BivariateType.PEARSONCORRELATION) {
                    if(variable0.equals(variable1)) {
                        value=1.0;
                    }
                    else {
                        value = PearsonCorrelation.calculateCorrelation(transposeDataList);
                    }
                }
                else if(type==BivariateType.SPEARMANCORRELATION) {
                    if(variable0.equals(variable1)) {
                        value=1.0;
                    }
                    else {
                        value = SpearmanCorrelation.calculateCorrelation(transposeDataList);
                    }
                }
                else if(type==BivariateType.KENDALLTAUCORRELATION) {
                    if(variable0.equals(variable1)) {
                        value=1.0;
                    }
                    else {
                        value = KendallTauCorrelation.calculateCorrelation(transposeDataList);
                    }
                }
                
                //bivariateMatrix.internalData.get(variable0).internalData.put(variable1, value);
                bivariateMatrix.put2d(variable0, variable1, value);
                
                if(!variable0.equals(variable1)) {
                    /*
                    if(!bivariateMatrix.internalData.containsKey(variable1)) {
                        bivariateMatrix.internalData.put(variable1, new AssociativeArray());
                    }
                    bivariateMatrix.internalData.get(variable1).internalData.put(variable0, value);
                    */
                    bivariateMatrix.put2d(variable1, variable0, value);
                }
            }
            //transposeDataList = null;
        }
        
        return bivariateMatrix;
    }
573691
private double variogram( double c0, double a, double sill, double rx, double ry, double rz ) {
        if (isNovalue(rz)) {
            rz = 0;
        }
        double value = 0;
        double h2 = Math.sqrt(rx * rx + rz * rz + ry * ry);
        if (pSemivariogramType == 0) {
            value = c0 + sill * (1 - Math.exp(-(h2 * h2) / (a * a)));
        }
        if (pSemivariogramType == 1) {
            // primotest semivariogram
            value = c0 + sill * (1 - Math.exp(-(h2) / (a)));
        }
        return value;
    }
53586
def _mean_square_error(y, y_pred, w):
    """"""
    return np.average(((y_pred - y) ** 2), weights=w)
669111
def MeanVarEstEQ(y, x, covariates, tol=1e-8):
    

    pcount = covariates.shape[0] + 2
    N = y.shape[0]
    beta_count = pcount * 2
    X = numpy.ones((N, pcount))
    X[:, 1] = x

    aprime = [1, 0]
    idx = 2
    for cov in covariates:
        aprime.append(-numpy.mean(cov)/numpy.std(cov))
        X[:, idx] = cov
        idx += 1
    aprime = numpy.matrix(aprime)

    # http://stackoverflow.com/questions/7267226/range-for-floats
    def frange(x, y, jump):
        while x < y:
            yield x
            x += jump
        yield 0

    class PhiReturn(object):
        def __init__(self, phi, dtheta):
            self.phi = phi
            self.dtheta = dtheta

    def dot_diag(x, y):
        """This is a conveniance function to perform the following \
        in a more efficient manner:

        :param x: arr1
        :param y: arr2

        x.dot(numpy.diag(y))

        y must be a single dimensioned array

        """
        if len(y.shape) != 1:
            print >> sys.stderr, "You can't pass an array of shape %s to dot_diag" % (y.shape)
            sys.exit(1)

        result = numpy.empty(x.shape)
        for i in range(0, x.shape[0]):
            result[i] = x[i] * y

        return result

    def diag_dot(x, y):
        """diagonal multiply
        """
        result = numpy.empty(y.shape)
        for i in range(0, y.shape[0]):
            result[i] = y[i] * x[i]
        return result

    def Phi(theta):
        MM = y - numpy.dot(X, theta[0:pcount])
        SS = numpy.exp(numpy.dot(-X, theta[pcount:]))
        DD1 = MM * SS
        DD2 = 0.5 * MM**2 * SS
        Phi = numpy.hstack((DD1.dot(X), (DD2-0.5).dot(X)))

        tX = X.transpose()

        t1 = numpy.empty(tX.shape)
        t2 = numpy.empty(tX.shape)
        t3 = numpy.empty(tX.shape)
        for i in range(0,tX.shape[0]):
            t1[i] = tX[i]*-SS
            t2[i] = tX[i]*-DD1
            t3[i] = tX[i]*-DD2
        t1 = t1.dot(X)
        t2 = t2.dot(X)
        t3 = t3.dot(X)

        dtheta = numpy.hstack((numpy.vstack((t1, t2)), numpy.vstack((t2, t3))))
        return PhiReturn(Phi, dtheta)



    class MvSolveReturn(object):
        def __init__(self, theta, dtheta):
            self.theta = theta
            self.dtheta = dtheta

    def MVsolve(theta_new):
        solution_found = False

        mvsolve_iterations = 0

        while not solution_found:
            theta_old = theta_new.copy()
            tmp = Phi(theta_old)
            #print "ITR", mvsolve_iterations, numpy.mean(theta_new),numpy.sum(numpy.absolute(theta_new - theta_old))
            theta_new = theta_old - scipy.linalg.solve(tmp.dtheta, tmp.phi)
            mvsolve_iterations += 1

            if (numpy.sum(numpy.absolute(theta_new - theta_old)) < tol):
                tmp = Phi(theta_new)
                solution_found = True
            else:
                if mvsolve_iterations > 25000:
                    #print >> sys.stderr, mvsolve_iterations, "failures"
                    raise UnsolvedLocus("")

        return MvSolveReturn(theta_new, tmp.dtheta), mvsolve_iterations

    def MVcalcB(theta):
        MM = y - X.dot(theta[0:pcount])
        SS = numpy.exp(-X.dot(theta[pcount:]))
        DD1 = MM * SS
        DD2 = 0.5 * MM**2 * SS
        AA = numpy.hstack((diag_dot(DD1, X), diag_dot(DD2-0.5,X)))

        return numpy.transpose(AA).dot(AA)/ N


    mod = None
    itr = 0
    total_iterations = 0
    for i in frange(0.00, 1.0, 0.05):
        theta=numpy.empty((beta_count))
        theta[:] = i

        try:
            mod, iterations = MVsolve(theta)
            total_iterations += iterations
            if i > 0.05:
                print >> sys.stderr, "Completed: ", total_iterations, itr
            break
        except exceptions.ValueError as e:
            pass
        except numpy.linalg.linalg.LinAlgError as e:
            pass
        except Exception as inst:
            #print type(inst)
            pass
        itr += 1

    if not mod:
        raise UnsolvedLocus("")

    try:
        ainv = scipy.linalg.inv(mod.dtheta) * N
    except:
        raise ValueError("Singular Matrix Encountered")
    B = MVcalcB(mod.theta)
    V = ainv.dot(B).dot(ainv.transpose())



    # Focus on the two parameters of interest
    theta2 = numpy.array([mod.theta[1], mod.theta[pcount+1]])
    V2 = V[1:beta_count:pcount,1:beta_count:pcount]
    pvalt = 1 - scipy.stats.chi2.cdf(theta2.dot(scipy.linalg.inv(V2)).dot(theta2) * N, 2)


    ## From Chun's updated code:
    theta= mod.theta
    se = numpy.sqrt(numpy.diag(V)/N)
    pval = 2*scipy.stats.norm.cdf(-numpy.absolute(mod.theta/se))
    return pvalt, theta, pval, se, V/N
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))
66037
def psicomputations(variance, Z, variational_posterior, return_psi2_n=False):
    
    # here are the "statistics" for psi0, psi1 and psi2
    # Produced intermediate results:
    # psi0    N
    # psi1    NxM
    # psi2    MxM
    mu = variational_posterior.mean
    S = variational_posterior.variance

    psi0 = (variance*(np.square(mu)+S)).sum(axis=1)
    Zv = variance * Z
    psi1 = np.dot(mu,Zv.T)
    if return_psi2_n:
        psi2 = psi1[:,:,None] * psi1[:,None,:] + np.dot(S[:,None,:] * Zv[None,:,:], Zv.T)
    else:
        psi2 = np.dot(S.sum(axis=0) * Zv, Zv.T) + tdot(psi1.T)

    return psi0, psi1, psi2
70764
function rf(df1, df2) {
      var rchisqRatio;

      rchisqRatio = function(d) { return chisq.rchisq(d)() / d; };

      if (df1 <= 0 || df2 <= 0) { return function() { return NaN; }; }

      return function() {
         return rchisqRatio(df1) / rchisqRatio(df2);
      };
   }
42562
def multi_reciprocal_extra(xs, ys, noise=False):
    
    ns = np.linspace(0.5, 6.0, num=56)
    best = ['', np.inf]
    fit_results = {}
    weights = get_weights(xs, ys)
    for n in ns:
        popt = extrapolate_reciprocal(xs, ys, n, noise)
        m = measure(reciprocal, xs, ys, popt, weights)
        pcov = []
        fit_results.update({n: {'measure': m, 'popt': popt, 'pcov': pcov}})
    for n in fit_results:
        if fit_results[n]['measure'] <= best[1]:
            best = reciprocal, fit_results[n]['measure'], n
    return fit_results[best[2]]['popt'], fit_results[best[2]]['pcov'], best
1023255
def ts_stats_significance(ts, ts_stat_func, null_ts_func, B=1000, permute_fast=False):
    
    stats_ts = ts_stat_func(ts)
    if permute_fast:
        # Permute it in 1 shot
        null_ts = map(np.random.permutation, np.array([ts, ] * B))
    else:
        null_ts = np.vstack([null_ts_func(ts) for i in np.arange(0, B)])
    stats_null_ts = np.vstack([ts_stat_func(nts) for nts in null_ts])
    pvals = []
    nums = []
    for i in np.arange(0, len(stats_ts)):
        num_samples = np.sum((stats_null_ts[:, i] >= stats_ts[i]))
        nums.append(num_samples)
        pval = num_samples / float(B)
        pvals.append(pval)

    return stats_ts, pvals, nums