• API
• FAQ
• Tools
• Trends
• Archive
daily pastebin goal
40%
SHARE
TWEET

# Untitled

a guest Sep 13th, 2017 43 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #%% Define model
2. import numpy,math
3. import matplotlib.pyplot as plt
4. import random as random
5.
6. import theano.tensor as t
7. import theano
8.
9. random.seed(1)  # set random seed
10.
11. # copy-pasted function from the specific model used in the source and adapted with as_op
12. @theano.compile.ops.as_op(itypes=[t.iscalar, t.dscalar, t.fscalar, t.fscalar],otypes=[t.dvector])
13. def sampleFromSalpeter(N, alpha, M_min, M_max):
14.     log_M_Min = math.log(M_min)
15.     log_M_Max = math.log(M_max)
16.     maxlik = math.pow(M_min, 1.0 - alpha)
17.     Masses = []
18.     while (len(Masses) < N):
19.         logM = random.uniform(log_M_Min,log_M_Max)
20.         M    = math.exp(logM)
21.         likelihood = math.pow(M, 1.0 - alpha)
22.         u = random.uniform(0.0,maxlik)
23.         if (u < likelihood):
24.             Masses.append(M)
25.     return Masses
26.
27. # SAME function as above, used to make test data (so no Theano here)
28. def sampleFromSalpeterDATA(N, alpha, M_min, M_max):
29.     log_M_Min = math.log(M_min)
30.     log_M_Max = math.log(M_max)
31.     maxlik = math.pow(M_min, 1.0 - alpha)
32.     Masses = []
33.     while (len(Masses) < N):
34.         logM = random.uniform(log_M_Min,log_M_Max)
35.         M    = math.exp(logM)
36.         likelihood = math.pow(M, 1.0 - alpha)
37.         u = random.uniform(0.0,maxlik)
38.         if (u < likelihood):
39.             Masses.append(M)
40.     return Masses
41.
42. # Generate toy data.
43. N      = 1000000  # Draw 1 Million stellar masses.
44. alpha  = 2.35
45. M_min  = 1.0
46. M_max  = 100.0
47. Masses = sampleFromSalpeterDATA(N, alpha, M_min, M_max)
48.
49. #%% Estimation process
50. import pymc3 as pm
51. basic_model = pm.Model()
52. with basic_model:
53.
54.     # Priors for unknown model parameters
55.     alpha2 = pm.Normal('alpha2', mu=3, sd=10)
56.
57.     N2=t.constant(1000000)
58.     M_min2  = t.constant(1.0)
59.     M_max2  = t.constant(100.0)
60.
61.     # Expected value of outcome
62.     m =  sampleFromSalpeter(N2, alpha2, M_min2, M_max2)
63.
64.     # Likelihood (sampling distribution) of observations
65.     Y_obs = pm.Normal('Y_obs', mu=m, sd=10, observed=Masses)
66.
67. #%% Sample
68. with basic_model:
69.     step = pm.Metropolis()
70.     trace = pm.sample(10000, step=step)
71.
72. #%% Plot posteriors
73. _ = pm.traceplot(trace)
74. pm.summary(trace)
RAW Paste Data
Top