mjam03

# define dist
my_dist = cust_dist(0.8)
# create samples
n = 1000000
cust_samps = my_dist.rvs(size=n)

# plot histogram
fig, ax = plt.subplots(ncols=1, figsize=(25,12))
# plot hist and pdf to see if they aline
sns.distplot(cust_samps, label='cust_samps', kde=False, norm_hist=True, ax=ax)

mjam03

# define our gaussian look-a-like distribution
class cust_dist(stats.rv_continuous):
    
    # define init with sigma deviation param e
    def __init__(self, e, *args, **kwargs):
        super().__init__(*args, **kwargs)
        # init our variance divergence
        self.e = e
        # init our cdf and ppf functions
        self.cdf_func, self.ppf_func = self.create_cdf_ppf()

mjam03

def _ppf_single(self, q, *args):
    factor = 10.
    left, right = self._get_support(*args)

    if np.isinf(left):
        left = min(-factor, right)
        while self._ppf_to_solve(left, q, *args) > 0.:
            left, right = left * factor, left
        # left is now such that cdf(left) <= q
        # if right has changed, then cdf(right) > q

mjam03

# define instance of dist
naive_dist = naive_cust_dist(e=0.5)
# draw samples
n = 1000
naive_samps = naive_dist.rvs(size=n)

# plot them
fig, ax = plt.subplots(ncols=1, figsize=(25,12))
# plot hist and pdf to see if they aline
sns.distplot(naive_samps, label='naive_samps', kde=False, norm_hist=True, ax=ax)

mjam03

# define our custom distribution
class naive_cust_dist(stats.rv_continuous):
    
    # define init with sigma deviation param e
    def __init__(self, e, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.e = e
    
    # function to return normal distribution pdf
    def norm_p(self, x, loc=0, scale=1):

mjam03

# usual suspects
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import scipy.stats as stats
import warnings

from scipy.integrate import simps
from scipy.interpolate import interp1d

mjam03

# sample size
n = 1000000
# define scipy implementation
snorm = stats.norm()
# define newer numpy implementation
nnorm = np.random.default_rng()

mjam03

double legacy_gauss(aug_bitgen_t *aug_state) {
  if (aug_state->has_gauss) {
    const double temp = aug_state->gauss;
    aug_state->has_gauss = false;
    aug_state->gauss = 0.0;
    return temp;
  } else {
    double f, x1, x2, r2;

    do {

mjam03

def __init__(self, seed=None):
    if seed is None:
        bit_generator = _MT19937()
    elif not hasattr(seed, 'capsule'):
        bit_generator = _MT19937()
        bit_generator._legacy_seeding(seed)
    else:
        bit_generator = seed

mjam03

# stripped out the doc string for ease
def check_random_state(seed):
    if seed is None or seed is np.random:
        return np.random.mtrand._rand
    
    ### this is the key line ###
    if isinstance(seed, (numbers.Integral, np.integer)):
        return np.random.RandomState(seed)
    
    if isinstance(seed, (np.random.RandomState, np.random.Generator)):