elksie5000
elksie5000

Reputation: 7752

How to implement simple Monte Carlo function in pymc

I'm trying to get my head around how to implement a Monte Carlo function in python using pymc to replicate a spreadsheet by Douglas Hubbard in his book How to Measure Anything

My attempt was:

import numpy as np
import pandas as pd
from pymc import DiscreteUniform, Exponential, deterministic, Poisson, Uniform, Normal, Stochastic, MCMC, Model


maintenance_saving_range = DiscreteUniform('maintenance_saving_range', lower=10, upper=21)
labour_saving_range = DiscreteUniform('labour_saving_range', lower=-2, upper=9)
raw_material_range = DiscreteUniform('maintenance_saving_range', lower=3, upper=10)
production_level_range = DiscreteUniform('maintenance_saving_range', lower=15000, upper=35000)



@deterministic(plot=False)
def rate(m = maintenance_saving_range, l = labour_saving_range, r=raw_material_range, p=production_level_range):
    return (m + l + r) * p

model = Model([rate, maintenance_saving_range, labour_saving_range, raw_material_range, production_level_range])
mc = MCMC(model)

Unfortunately, I'm getting an error: ValueError: A tallyable PyMC object called maintenance_saving_range already exists. This will cause problems for some database backends.

What have I got wrong?

Upvotes: 2

Views: 290

Answers (1)

elksie5000
elksie5000

Reputation: 7752

Ah, it was a copy and paste error.

I'd called three distributions by the same name.

Here's the code that works.

import numpy as np
import pandas as pd
from pymc import DiscreteUniform, Exponential, deterministic, Poisson, Uniform, Normal, Stochastic, MCMC, Model
%matplotlib inline
import matplotlib.pyplot as plt

maintenance_saving_range = DiscreteUniform('maintenance_saving_range', lower=10, upper=21)
labour_saving_range = DiscreteUniform('labour_saving_range', lower=-2, upper=9)
raw_material_range = DiscreteUniform('raw_material_range', lower=3, upper=10)
production_level_range = DiscreteUniform('production_level_range', lower=15000, upper=35000)


@deterministic(plot=False, name="rate")
def rate(m = maintenance_saving_range, l = labour_saving_range, r=raw_material_range, p=production_level_range):
    #out = np.empty(10000)
    out = (m + l + r) * p
    return out
model = Model([rate, maintenance_saving_range, labour_saving_range, raw_material_range])
mc = MCMC(model)
mc.sample(iter=10000)

Upvotes: 1

Related Questions