def pca(X, nb_components=0):

    mu=np.mean(X,axis=0)
    #subtract mean from data
    Y=X.copy()
    for row in range(0,28):
        Y[row,:]-=mu

    #calculate covariance matrix
    C=np.cov(Y,rowvar=False)

    Eigva,Eigve=np.linalg.eigh(C)
    sort_indices = Eigva.argsort()[::-1]                             # getting their correct order - decreasing

    Eigva = Eigva[sort_indices]                               # puttin the evalues in that order
    Eigve = Eigve[:,sort_indices]                             # same for the evectors

    totvar1=np.sum(Eigva)

    Eigva=Eigva[0:nb_components]
    Eigve=Eigve[:,0:nb_components]
    print "eigva"
    print Eigva
    print
    totvar2=np.sum(Eigva)
    verh=totvar2/totvar1
    print verh

    norms = np.linalg.norm(Eigve, axis=0)                           # find the norm of each eigenvector

    Eigve = Eigve / norms                                   # normalize all eigenvectors

    b=np.dot(Eigve.T,Y.T)
    print b
    X2=np.dot(b.T,Eigve.T)

    for row in range(0,28):
        X2[row,:]+=mu

    #de B-parameters aanpassen en testen op het mean model
    b2=np.zeros([nb_components,1])
    b2[2]=0+2*np.sqrt(Eigva[2])
    mu2=np.dot(Eigve,b2)
    mu2=mu2.flatten()
    mu2+=mu
    print mu2
    print mu2.shape
    return X2,mu2