Выборка Гиббса с полиномиальной моделью не может обновить параметр - PullRequest
Новая
       15

Выборка Гиббса с полиномиальной моделью не может обновить параметр

0 голосов
/ 24 января 2019

Я пытаюсь реализовать алгоритм, взятый у Альберта и Чиба, для выборки со скрытым полиномиальным распределением.Кажется, код работает до тех пор, пока я не попытаюсь обновить ковариационную матрицу: тогда я не могу сгенерировать нормальное случайное число. 

def gibbs_u(x,y,sim):

    x = x.values #need transform as it is a dataframe
    y = y.values #same

    n,k = x.shape


    #Get true number of categories:
    J = np.max(y)+1

    #initialize beta to MLE values and rho to 0:
    rho = 0.01
    beta = 0.238

    #initialize Z to empty array
    Z = np.empty((n,J))

    #conditional distribution of Z:
    cov_mat = np.identity(J)
    cov_mat[0,1] = rho
    cov_mat[1,0] = rho

    betaSim = np.empty((sim+1, 1))  #Matrice de retour avec tous les beta simulés
    betaSim[0,] = beta
    rhoSim = np.empty((sim+1, 1))  #Matrice de retour avec tous les rho simulés
    rhoSim[0,] = rho


    for i in range(sim):
        #Update Z
        for j in range(n):
            x = np.reshape(x, (n,J), order='C')
            Z = np.reshape(Z, (n,J), order='C')
            meanZ = x[j,:]*beta
            while True:
                Z[j,:] = np.random.multivariate_normal(meanZ.T, cov_mat)
                if np.argmax(Z[j,:]) == y[j]:
                    break

        # Update beta:
        x = np.reshape(x, (n*J, 1), order='C')
        Z = np.reshape(Z, (n*J, 1), order='C')
        Omega = np.kron(np.identity(n), cov_mat)
        xTOx = np.linalg.inv(np.dot(x.T,np.dot(np.linalg.inv(Omega), x)))
        beta_Z = np.dot(xTOx,np.dot(x.T, np.dot(np.linalg.inv(Omega), Z)))
        beta = np.random.normal(beta_Z[0][0], xTOx)[0][0]

        #Update rho:
        w = Z - x*beta
        S_w1 = np.sum((w[::3]-np.mean(w[::3]))**2)
        S_w1w2 = np.sum((w[::3]-np.mean(w[::3]))*(w[1::3]-np.mean(w[1::3])))
        S_w2 = np.sum((w[1::3]-np.mean(w[1::3]))**2)
        #from Lee 1989, chap.6:
        r = S_w1w2/(np.sqrt(S_w1*S_w2))
        rho = (1-rho**2)**((n-1)/2)/(1-rho*r)**(n-(3/2))
        print(rho)

        # and of covariance w/ new rho: 
        #cov_mat[0,1] = round(rho, 4)
        #cov_mat[1,0] = round(rho, 4)

        betaSim[i+1,] = beta
        rhoSim[i+1,] = rho


    return(betaSim, rhoSim)

Я попытался использовать округление rho, чтобы убедиться, что ковариация пригодна для использования, ноэто не похоже на работу.Видите ли вы другую стратегию, которую я мог бы использовать?

...