pymc3:具有多维浓度因子的狄利克雷

我正在努力实现一个模型,其中狄利克雷变量的集中因子取决于另一个变量。


情况如下:


系统因组件故障而失败(有三个组件,每次测试/观察只有一个失败)。


组件失效的概率取决于温度。


这是该情况的(评论)简短实现:


import numpy as np

import pymc3 as pm

import theano.tensor as tt



# Temperature data : 3 cold temperatures and 3 warm temperatures

T_data = np.array([10, 12, 14, 80, 90, 95])


# Data of failures of 3 components : [0,0,1] means component 3 failed

F_data = np.array([[0, 0, 1], \

       [0, 0, 1], \

       [0, 0, 1], \

       [1, 0, 0], \

       [1, 0, 0], \

       [1, 0, 0]])


n_component = 3


# When temperature is cold : Component 1 fails

# When temperature is warm : Component 3 fails

# Component 2 never fails


# Number of observations :

n_obs = len(F_data)



# The number of failures can be modeled as a Multinomial F ~ M(n_obs, p) with parameters 

# -  n_test : number of tests (Fixed)

# -  p : probability of failure of each component (shape (n_obs, 3))


# The probability of failure of components follows a Dirichlet distribution p ~ Dir(alpha) with parameters:

# -  alpha : concentration (shape (n_obs, 3))

# The Dirichlet distributions ensures the probabilities sum to 1 


# The alpha parameters (and the the probability of failures) depend on the temperature alpha ~ a + b * T

# - a : bias term (shape (1,3))

# - b : describes temperature dependency of alpha (shape (1,3))

_


# The prior on "a" is a normal distributions with mean 1/2 and std 0.001

# a ~ N(1/2, 0.001)


# The prior on "b" is a normal distribution zith mean 0 and std 0.001

# b ~ N(0, 0.001)



# Coding it all with pymc3


with pm.Model() as model:

    a = pm.Normal('a', 1/2, 1/(0.001**2), shape = n_component)

    b = pm.Normal('b', 0, 1/(0.001**2), shape = n_component)


    # I generate 3 alphas values (corresponding to the 3 components) for each of the 6 temperatures

    # I tried different ways to compute alpha but nothing worked out

    alphas = pm.Deterministic('alphas', a + b * tt.stack([T_data, T_data, T_data], axis=1))

    #alphas = pm.Deterministic('alphas', a + b[None, :] * T_data[:, None])

    #alphas = pm.Deterministic('alphas', a + tt.outer(T_data,b))




慕少森
浏览 132回答 1
1回答
打开App,查看更多内容
随时随地看视频慕课网APP

相关分类

Python