因前一篇https://blog.csdn.net/fjssharpsword/article/details/97000479采樣問題未解決，發(fā)現(xiàn)如下github上有BPMF代碼，采用wishart先驗，性能和pymc3一致。

參考：https://github.com/LoryPack/BPMF

# coding:utf-8 ''' @author: Jason.F @data: 2019.08.01 @function: baseline BPMF(Bayesian Probabilistic Matrix Factorization)Datatset: MovieLens-1m:https://grouplens.org/datasets/movielens/ Evaluation: RMSE '''import numpy as np import random import pandas as pd from numpy.random import multivariate_normal from scipy.stats import wishartclass DataSet:def __init__(self):self.trainset, self.testset, self.maxu, self.maxi, self.maxr = self._getDataset_as_list()def _getDataset_as_list(self):#trainsetfilePath = "/data/fjsdata/BMF/ml-1m.train.rating" data = pd.read_csv(filePath, sep='\t', header=None, names=['user', 'item', 'rating'], \usecols=[0, 1, 2], dtype={0: np.int32, 1: np.int32, 2: np.float})maxu, maxi, maxr = data['user'].max()+1, data['item'].max()+1, data['rating'].max()print('Dataset Statistics: Interaction = %d, User = %d, Item = %d, Sparsity = %.4f' % \(data.shape[0], maxu, maxi, data.shape[0]/(maxu*maxi)))trainset = data.values.tolist()#testsetfilePath = "/data/fjsdata/BMF/ml-1m.test.rating" data = pd.read_csv(filePath, sep='\t', header=None, names=['user', 'item', 'rating'], \usecols=[0, 1, 2], dtype={0: np.int32, 1: np.int32, 2: np.float})testset = data.values.tolist()return trainset, testset, maxu, maxi, maxr def list_to_matrix(self, dataset, maxu, maxi): dataMat = np.zeros([maxu, maxi], dtype=np.float32)for u,i,r in dataset:dataMat[int(u)][int(i)] = float(r)return np.array(dataMat)def Normal_Wishart(mu_0, lamb, W, nu, seed=None):"""Function extracting a Normal_Wishart random variable"""# first draw a Wishart distribution:Lambda = wishart(df=nu, scale=W, seed=seed).rvs() # NB: Lambda is a matrix.# then draw a Gaussian multivariate RV with mean mu_0 and(lambda*Lambda)^{-1} as covariance matrix.cov = np.linalg.inv(lamb * Lambda) # this is the bottleneck!!mu = multivariate_normal(mu_0, cov)return mu, Lambda, covdef BPMF(R, R_test, U_in, V_in, T, D, initial_cutoff, lowest_rating, highest_rating,mu_0=None, Beta_0=None, W_0=None, nu_0=None):"""R is the ranking matrix (NxM, N=#users, M=#movies); we are assuming that R[i,j]=0 means that user i has not ranked movie jR_test is the ranking matrix that contains test values. Same assumption as above. U_in, V_in are the initial values for the MCMC procedure. T is the number of steps. D is the number of hidden features that are assumed in the model. mu_0 is the average vector used in sampling the multivariate normal variableBeta_0 is a coefficient (?)W_0 is the DxD scale matrix in the Wishart sampling nu_0 is the number of degrees of freedom used in the Wishart sampling. U matrices are DxN, while V matrices are DxM."""def ranked(i, j): # function telling if user i ranked movie j in the train dataset.if R[i, j] != 0:return Trueelse:return Falsedef ranked_test(i, j): # function telling if user i ranked movie j in the test dataset.if R_test[i, j] != 0:return Trueelse:return FalseN = R.shape[0]M = R.shape[1]R_predict = np.zeros((N, M))U_old = np.array(U_in)V_old = np.array(V_in)train_err_list = []test_err_list = []train_epoch_list = []# initialize now the hierarchical priors:alpha = 2 # observation noise, they put it = 2 in the papermu_u = np.zeros((D, 1))mu_v = np.zeros((D, 1))Lambda_U = np.eye(D)Lambda_V = np.eye(D)# COUNT HOW MAY PAIRS ARE IN THE TEST AND TRAIN SET:pairs_test = 0pairs_train = 0for i in range(N):for j in range(M):if ranked(i, j):pairs_train = pairs_train + 1if ranked_test(i, j):pairs_test = pairs_test + 1# print(pairs_test, pairs_train)# SET THE DEFAULT VALUES for Wishart distribution# we assume that parameters for both U and V are the same.if mu_0 is None:mu_0 = np.zeros(D)if nu_0 is None:nu_0 = Dif Beta_0 is None:Beta_0 = 2if W_0 is None:W_0 = np.eye(D)# results = pd.DataFrame(columns=['step', 'train_err', 'test_err'])for t in range(T):# print("Step ", t)# FIRST SAMPLE THE HYPERPARAMETERS, conditioned on the present step user and movie feature matrices U_t and V_t:# parameters common to both distributions:Beta_0_star = Beta_0 + Nnu_0_star = nu_0 + NW_0_inv = np.linalg.inv(W_0) # compute the inverse once and for all# movie hyperparameters:V_average = np.sum(V_old, axis=1) / N # in this way it is a 1d array!!# print (V_average.shape)S_bar_V = np.dot(V_old, np.transpose(V_old)) / N # CHECK IF THIS IS RIGHT!mu_0_star_V = (Beta_0 * mu_0 + N * V_average) / (Beta_0 + N)W_0_star_V_inv = W_0_inv + N * S_bar_V + Beta_0 * N / (Beta_0 + N) * np.dot(np.transpose(np.array(mu_0 - V_average, ndmin=2)), np.array((mu_0 - V_average), ndmin=2))W_0_star_V = np.linalg.inv(W_0_star_V_inv)mu_V, Lambda_V, cov_V = Normal_Wishart(mu_0_star_V, Beta_0_star, W_0_star_V, nu_0_star, seed=None)# user hyperparameters# U_average=np.transpose(np.array(np.sum(U_old, axis=1)/N, ndmin=2)) #the np.array and np.transpose are needed for it to be a column vectorU_average = np.sum(U_old, axis=1) / N # in this way it is a 1d array!! #D-long# print (U_average.shape)S_bar_U = np.dot(U_old, np.transpose(U_old)) / N # CHECK IF THIS IS RIGHT! #it is DxDmu_0_star_U = (Beta_0 * mu_0 + N * U_average) / (Beta_0 + N)W_0_star_U_inv = W_0_inv + N * S_bar_U + Beta_0 * N / (Beta_0 + N) * np.dot(np.transpose(np.array(mu_0 - U_average, ndmin=2)), np.array((mu_0 - U_average), ndmin=2))W_0_star_U = np.linalg.inv(W_0_star_U_inv)mu_U, Lambda_U, cov_U = Normal_Wishart(mu_0_star_U, Beta_0_star, W_0_star_U, nu_0_star, seed=None)# print (S_bar_U.shape, S_bar_V.shape)# print (np.dot(np.transpose(np.array(mu_0-U_average, ndmin=2)),np.array((mu_0-U_average), ndmin=2).shape))# UP TO HERE IT PROBABLY WORKS, FROM HERE ON IT HAS TO BE CHECKED!!!"""SAMPLE THEN USER FEATURES (possibly in parallel):"""U_new = np.array([]) # define the new stuff.V_new = np.array([])for i in range(N): # loop over the users# first compute the parameters of the distributionLambda_U_2 = np.zeros((D, D)) # second term in the construction of Lambda_Umu_i_star_1 = np.zeros(D) # first piece of mu_i_starfor j in range(M): # loop over the moviesif ranked(i, j): # only if movie j has been ranked by user i!Lambda_U_2 = Lambda_U_2 + np.dot(np.transpose(np.array(V_old[:, j], ndmin=2)),np.array((V_old[:, j]), ndmin=2)) # CHECKmu_i_star_1 = V_old[:, j] * R[i, j] + mu_i_star_1 # CHECK DIMENSIONALITY!!!!!!!!!!!!# coeff=np.transpose(np.array(V_old[j]*R[i,j], ndmin=2))+coeff #CHECK DIMENSIONALITY!!!!!!!!!!!!Lambda_i_star_U = Lambda_U + alpha * Lambda_U_2Lambda_i_star_U_inv = np.linalg.inv(Lambda_i_star_U)mu_i_star_part = alpha * mu_i_star_1 + np.dot(Lambda_U,mu_U) ###CAREFUL!! Multiplication matrix times a row vector!! It should give as an output a row vector as for how it worksmu_i_star = np.dot(Lambda_i_star_U_inv, mu_i_star_part)# extract now the U values!U_new = np.append(U_new, multivariate_normal(mu_i_star, Lambda_i_star_U_inv))# you need to reshape U_new and transpose it!!U_new = np.transpose(np.reshape(U_new, (N, D)))# print (U_new.shape)"""SAMPLE THEN MOVIE FEATURES (possibly in parallel):"""for j in range(M):Lambda_V_2 = np.zeros((D, D)) # second term in the construction of Lambda_Umu_i_star_1 = np.zeros(D) # first piece of mu_i_starfor i in range(N): # loop over the moviesif ranked(i, j):Lambda_V_2 = Lambda_V_2 + np.dot(np.transpose(np.array(U_new[:, i], ndmin=2)),np.array((U_new[:, i]), ndmin=2))mu_i_star_1 = U_new[:, i] * R[i, j] + mu_i_star_1 # CHECK DIMENSIONALITY!!!!!!!!!!!!# coeff=np.transpose(np.array(V_old[j]*R[i,j], ndmin=2))+coeff #CHECK DIMENSIONALITY!!!!!!!!!!!!Lambda_j_star_V = Lambda_V + alpha * Lambda_V_2Lambda_j_star_V_inv = np.linalg.inv(Lambda_j_star_V)mu_i_star_part = alpha * mu_i_star_1 + np.dot(Lambda_V, mu_V)mu_j_star = np.dot(Lambda_j_star_V_inv, mu_i_star_part)V_new = np.append(V_new, multivariate_normal(mu_j_star, Lambda_j_star_V_inv))# you need to reshape U_new and transpose it!!V_new = np.transpose(np.reshape(V_new, (M, D)))# save U_new and V_new in U_old and V_old for next iteration: U_old = np.array(U_new)V_old = np.array(V_new)if t > initial_cutoff: # initial_cutoff is needed to discard the initial transientR_step = np.dot(np.transpose(U_new), V_new)for i in range(N): # reduce all the predictions to the correct ratings range.for j in range(M):if R_step[i, j] > highest_rating:R_step[i, j] = highest_ratingelif R_step[i, j] < lowest_rating:R_step[i, j] = lowest_ratingR_predict = (R_predict * (t - initial_cutoff - 1) + R_step) / (t - initial_cutoff)train_err = 0 # initialize the errors.test_err = 0# compute now the RMSE on the train dataset:for i in range(N):for j in range(M):if ranked(i, j):train_err = train_err + (R_predict[i, j] - R[i, j]) ** 2train_err_list.append(np.sqrt(train_err / pairs_train))print("Training RMSE at iteration ", t - initial_cutoff, " : ", "{:.4}".format(train_err_list[-1]))# compute now the RMSE on the test dataset:for i in range(N):for j in range(M):if ranked_test(i, j):test_err = test_err + (R_predict[i, j] - R_test[i, j]) ** 2test_err_list.append(np.sqrt(test_err / pairs_test))print("Test RMSE at iteration ", t - initial_cutoff, " : ", "{:.4}".format(test_err_list[-1]))return R_predictif __name__ == "__main__":ds = DataSet()#loading dataset\R = ds.list_to_matrix(ds.trainset, ds.maxu, ds.maxi)#get matrixR_test = ds.list_to_matrix(ds.testset, ds.maxu, ds.maxi)#get matrixfor K in [8, 16, 32, 64]:U_in = np.zeros((K, ds.maxu)) V_in = np.zeros((K, ds.maxi))R_pred = BPMF(R, R_test, U_in, V_in, T=100, D=K, initial_cutoff=0, lowest_rating=0, highest_rating=ds.maxr)squaredError = []for u,i,r in ds.testset:error=r - nR[int(u)][int(i)]squaredError.append(error * error)rmse =math.sqrt(sum(squaredError) / len(squaredError))print("RMSE@{}:{}".format(K, rmse))

?

總結(jié)

以上是生活随笔為你收集整理的推荐经典算法实现之BPMF(python+MovieLen)的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。