Построение процесса разрыва палки в R на основе кода Python - PullRequest
Новая
       73

Построение процесса разрыва палки в R на основе кода Python

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

Я хотел бы воспроизвести код Python в R-код о процессе взлома Stick, который является одной из схем построения для процесса Дирихле. Тем не менее, график, который я нарисовал внутри R, довольно отличается тем, что выборочные распределения DP не находятся вокруг базового распределения, H.

Ссылочный код Python взят из блога Остина Рочфорда . 

from matplotlib import pyplot as plt
import numpy as np
import pymc3 as pm
import scipy.stats as ss
import seaborn as sns
from statsmodels.datasets import get_rdataset
from theano import tensor as T
np.random.seed(433)
N=20    
K=30    
alpha=50
H = ss.norm # base dist
beta = ss.beta.rvs(1,alpha, size=(N,K))                          
pi = np.empty_like(beta)
pi[:, 0] = beta[:,0]
pi[:, 1:] = beta[:, 1:] * (1-beta[:, :-1]).cumprod(axis=1)
omega = H.rvs(size=(N,K))

x_plot = np.linspace(-3,3,200)            
sample_cdfs = (pi[..., np.newaxis]* np.less.outer(omega, x_plot)).sum(axis=1)

fig, ax = plt.subplots(figsize=(8,6))

ax.plot (x_plot, sample_cdfs[0],c="gray", alpha=0.75, label = "DP sample CDFs")
ax.plot(x_plot, sample_cdfs[1:].T, c="gray", alpha=0.75) 
ax.plot(x_plot, H.cdf(x_plot), c= "k", label = "Base CDF")

ax.set_title(r'$\alpha = {}$'.format(alpha))
ax.legend(loc=2)

Цифра справа является результатом в коде Python.

DP with alpha=50

И я попытался преобразовать его в код R: 

library(yarrr)

N=20;K=30;ngrid=200;alpha=50
xgrid = seq(-3,3,length.out=ngrid)
betas = matrix(rbeta(N*K, 1, alpha),nr=N, nc=K)
stick.to.right = c(1, cumprod(1 - betas))[1:K]
pis.temp = stick.to.right * betas
omega = matrix(rnorm(N*K),nr=N,nc=K)

dirac = array(numeric(N*K*ngrid),dim=c(N,K,ngrid))

for(i in 1:N){
    for(j in 1:K){
        for(k in 1:ngrid){
            dirac[i,j,k]=ifelse(omega[i,j]<xgrid[k],TRUE,FALSE)
        }
    }
}

pis = array(pis.temp,dim=c(N,K,200))

sample_cdfs = apply(pis* dirac,c(1,3),sum)

plot(xgrid,sample_cdfs[1,],col=piratepal("pony"),type="l",lwd=1,ylim=c(0,1))
for(i in 2:N) lines(xgrid,sample_cdfs[i,],col=piratepal("pony")[i])
lines(xgrid,pnorm(xgrid),lwd=2)

Я нарисовал сюжет DP с альфа = 50:

DP with alpha=50

Как я могу изменить код R, чтобы получить результат, аналогичный коду Python?

...