3

私は、いくつかの精神物理学的行動データの階層モデルを取得して、pymc3 で実行することに取り組んできました。私は全体的に信じられないほど感銘を受けましたが、Theano と pymc3 の速度を上げようとした後、ほとんど機能するモデルを手に入れましたが、いくつかの問題があります。

このコードは、パラメーター化されたバージョンのワイブルを 7 セットのデータに適合させるように作成されています。各試行は、バイナリ ベルヌーイの結果としてモデル化されますが、しきい値 (高さ、幅、高さ (典型的なガウスの a、c、および d) のガウス関数に適合するために使用される y 値としての出力)。

パラメータ化されたワイブルの使用はうまく機能しているようで、ワイブルの勾配が階層化され、データのチャンクごとにしきい値が個別に適合されます。ただし、 k および y_est から得られる出力は、それらが正しいサイズではない可能性があると思わせ、確率分布とは異なり、形状を指定できるようには見えません (これを行うための theano の方法がない限り)。私は見つけていません-しかし、私が読んだことから、theanoで形状を指定するのは難しいです)。

最終的には、y_est を使用してガウスの高さまたは幅を推定したいと思いますが、現在の出力は、y_est と k のサイズの問題に起因すると思われる信じられないほどの混乱をもたらします。以下のコードは、いくつかのデータをシミュレートする必要があり、その後にモデルが続きます。このモデルは、個々のしきい値をフィッティングして勾配を取得するのに優れた仕事をしますが、残りの部分を処理するとうまくいきません。

ご覧いただきありがとうございます。これまでのところ、pymc3 には非常に感銘を受けています。

編集：さて、y_est.tag.test_value.shape による形状出力は次のようになります

y_est.tag.test_value.shape
(101, 7)
k.tag.test_value.shape
(7,)

ここで問題が発生していると思いますが、私の側で構築が不十分なだけかもしれません。k は適切な形をしています (unique_xval ごとに 1 つの k 値)。y_est は、難易度ごとに 1 つの推定値 (unique_xval ごとに 1 つの y_est) ではなく、データ セット全体 (101x7) を出力しています。これを制御するために y_est が df_y_vals の特定のサブセットを取得するように指定する方法はありますか?

#Import necessary modules and define our weibull function
import numpy as np
import pylab as pl    
from scipy.stats import bernoulli

#x stimulus intensity
#g chance (0.5 for 2AFC)
# m slope
# t threshold
# a performance level defining threshold 
def weib(x,g,a,m,t):
    k=-np.log(((1-a)/(1-g))**(1/t))
    return 1- (1-g)*np.exp(- (k*x/t)**m);

#Output values from weibull function
xit=101
xvals=np.linspace(0.05,1,xit)
out_weib=weib(xvals, 0.5, 0.8, 3, 0.6)

#Okay, fitting the perfect output of a Weibull should be easy, contaminate         with some noise
#Slope of 3, threshold of 0.6


#How about 5% noise!

noise=0.05*np.random.randn(np.size(out_weib))
out=out_weib+noise

#Let's make this more like a typical experiment - 
#i.e. no percent correct, just one or zero
#Randomly pick based on the probability at each point whether they got the trial right or wrong
trial=np.zeros_like(out)
for i in np.arange(out.size):
    p=out_weib[i]
    trial[i] = bernoulli.rvs(p)

#Iterate for 6 sets of data, similar slope (from a normal dist), different thresh (output from gaussian)
#Gauss parameters=

true_gauss_height = 0.3
true_gauss_width = 0.01
true_gauss_elevation = 0.2

#What thresholds will we get then? 6 discrete points along that gaussian, from 0 to 180 degree mask

x_points=[0, 30, 60, 90, 120, 150, 180]

x_points=np.asarray(x_points)
gauss_points=true_gauss_height*np.exp(-    ((x_points**2)/2*true_gauss_width**2))+true_gauss_elevation

import pymc as pm2
import pymc3 as pm
import pandas as pd

slopes=pm2.rnormal(3, 3, size=7)
out_weib=np.zeros([xvals.size,x_points.size])

for i in np.arange(x_points.size):
    out_weib[:,i]=weib(xvals, 0.5, 0.8, slopes[i], gauss_points[i])

#Let's make this more like a typical experiment - i.e. no percent correct, just one or zero
#Randomly pick based on the probability at each point whether they got the trial right or wrong
trials=np.zeros_like(out_weib)

for i in np.arange(len(trials)):
    for ii in np.arange(gauss_points.size):
        p=out_weib[i,ii]
        trials[i,ii] = bernoulli.rvs(p)

#Let's make that data into a DataFrame for pymc3
y_vals=np.tile(xvals, [7, 1])

df_correct = pd.DataFrame(trials, columns=x_points)
df_y_vals = pd.DataFrame(y_vals.T, columns=x_points)
unique_xvals=x_points

import theano as th

with pm.Model() as hierarchical_model:
    # Hyperpriors for group node
    mu_slope = pm.Normal('mu_slope', mu=3, sd=1)
    sigma_slope = pm.Uniform('sigma_slope', lower=0.1, upper=2)

#Priors for the overall gaussian function - 3 params, the height of the gaussian
#Width, and elevation

gauss_width = pm.HalfNormal('gauss_width', sd=1)
gauss_elevation = pm.HalfNormal('gauss_elevation', sd=1)

slope = pm.Normal('slope', mu=mu_slope, sd=sigma_slope,     shape=unique_xvals.size)

thresh=pm.Uniform('thresh', upper=1, lower=0.1, shape=unique_xvals.size)

k = -th.tensor.log(((1-0.8)/(1-0.5))**(1/thresh))
y_est=1-(1-0.5)*th.tensor.exp(-(k*df_y_vals/thresh)**slope)

#We want our model to predict either height or width...height would be easier.
#Our Gaussian function has y values estimated by y_est as the 82% thresholds
#and Xvals based on where each of those psychometrics were taken.
#height_est=pm.Deterministic('height_est', (y_est/(th.tensor.exp((-unique_xvals**2)/2*gauss_width)))+gauss_elevation)
height_est = pm.Deterministic('height_est', (y_est-gauss_elevation)*th.tensor.exp((unique_xvals**2)/2*gauss_width**2))

#Define likelihood as Bernoulli for each binary trial
likelihood = pm.Bernoulli('likelihood',p=y_est, shape=unique_xvals.size, observed=df_correct)

#Find start
start=pm.find_MAP()
step=pm.NUTS(state=start)
#Do MCMC
trace = pm.sample(5000, step, njobs=1, progressbar=True) # draw 5000 posterior samples using NUTS sampling
4

1 に答える 1

0

「これを制御するために y_est が df_y_vals の特定のサブセットを取得するように指定する方法はありますか」と言ったときに、何をしたいのか正確にはわかりません。各 y_est 値について、使用することになっている df_y_vals の値を説明できますか? df_y_vals の形状は? y_est の形状はどのようなものですか? （7、）？

あなたが望むのは、 numpy の高度なインデックス作成を使用して df_y_vals にインデックスを作成することだと思います。これは、numpy と同様に PyMC でも機能します。詳細情報なしで正確に言うのは難しいです。

于 2015-05-08T18:35:08.493 に答える