How to optimize a model to joint multiple distributions

I’m trying to build a kind of volatility model.

import numpy as np
from pymc3.distributions.timeseries import GaussianRandomWalk
from scipy import stats
from scipy import optimize
import pymc3 as pm

x_mean = 10 
x_normal_dev = 10
x_outlier_dev = 30
x_normal_ratio = 0.70

binom = np.random.binomial(1, 0.7, 10000)
a = np.random.normal(x_mean, x_normal_dev, 10000)
b = np.random.normal(x_mean, x_outlier_dev, 10000)
c = a * binom + b * (1 - binom)

I made a sample data with the preceding code. Then, I wanted to make a model and estimate the normal_dev for the model.

with pm.Model() as model:
    #mu = pm.Uniform('mean', lower=-100, upper=100)
    mu = 10
    normal_dev = pm.Uniform('normal_dev', lower=0, upper=30)
    dev1 = pm.Normal('dev1', mu=mu, sigma=normal_dev)
    dev2 = pm.Normal('dev2', mu=mu, sigma=x_outlier_dev)
    binom = pm.Binomial('binom', 1, 0.7)
    output = pm.Normal('output', mu=dev1*binom + dev2*(1-binom),  tau=1,  observed=c)

with model:
    trace = pm.sample(1000)


It seems not to work well to estimate even only normal_dev value. What am I wrong with it?


It seems apply estimation not only normal_dev but others like binom value. Is there anyway to fix the variable value? I tried to use pm.Deterministic but it didn’t help.