SHARE
TWEET

Untitled

a guest Sep 13th, 2017 41 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