Python scipy.stats.norm.rvs() Examples

def weiner_process_fn(num_points, delta, x0=0, dt=1):
    """
    Generates Weiner process realisation trajectory.

    Args:
        num_points:     int, trajectory length;
        delta:          float, speed parameter;
        x0:             float, starting point;
        dt:             int, time increment;

    Returns:
        generated data as 1D np.array
    """
    x0 = np.asarray(x0)
    r = norm.rvs(size=x0.shape + (num_points,), scale=delta * (dt**.5))

    return np.cumsum(r, axis=-1) + np.expand_dims(x0, axis=-1) 

def impute(self, X):
        """Generate imputations using predictions from the fit linear model.

        The impute method returns the values for imputation. Missing values
        in a given dataset are replaced with the predictions from the least
        squares regression line of best fit plus a random draw from the normal
        error distribution.

        Args:
            X (pd.DataFrame): predictors to determine imputed values.

        Returns:
            np.array: imputed dataset.
        """
        # check if fitted then predict with least squares
        check_is_fitted(self, "statistics_")
        mse = self.statistics_["param"]
        preds = self.lm.predict(X)

        # add random draw from normal dist w/ mean squared error
        # from observed model. This makes lm stochastic
        mse_dist = norm.rvs(loc=0, scale=sqrt(mse), size=len(preds))
        imp = preds + mse_dist
        return imp 

def impute(self, feature=None, strategy=None, distribution=None):
        if self._matrix is None:
            raise ValueError('Must call set_input_matrix() before impute().')
        # If user does not specify which feature to impute, impute values
        # for all columns in the matrix.
        if feature is None:
            feature = FeatureMatrixTransform.ALL_FEATURES

        # If an imputation strategy is not specified, default to mean.
        if strategy is None:
            strategy = FeatureMatrixTransform.IMPUTE_STRATEGY_MEAN

        # If distribution is not specified, default to norm.
        if distribution is None:
            distribution = norm.rvs

        # TODO sxu: modify other modules to also return stuff
        if feature == FeatureMatrixTransform.ALL_FEATURES:
            self._impute_all_features(strategy, distribution=distribution)
        else:
            return self._impute_single_feature(feature, strategy, distribution=distribution) 

def _generatePos(self, lenBackground, lenSubstring, additionalInfo):
        from scipy.stats import norm
        center = (lenBackground-lenSubstring)/2.0
        validPos = False
        totalTries = 0
        while (validPos == False):
            sampledPos = int(norm.rvs(loc=center+self.offsetFromCenter,
                          scale=self.stdInBp))
            totalTries += 1
            if (sampledPos > 0 and sampledPos < (lenBackground-lenSubstring)):
                validPos = True
            if (totalTries%10 == 0 and totalTries > 0):
                print("Warning: made "+str(totalTries)+" attempts at sampling"
                      +" a position with lenBackground "+str(lenBackground)
                      +" and center "+str(center)+" and offset "
                      +str(self.offsetFromCenter)) 
        return sampledPos 

def test_aic_posterior_dependence():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    p2 = norm.logpdf(d, scale=2)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    c.add_chain(d, posterior=p2, num_free_params=1, num_eff_data_points=1000)
    aics = c.comparison.aic()
    assert len(aics) == 2
    assert aics[0] == 0
    expected = 2 * np.log(2)
    assert np.isclose(aics[1], expected, atol=1e-3) 

def setup(self):

        #########
        # PART 1: Make model calcium data
        #########

        # Data parameters
        RATE = 1       # mean firing rate of poisson spike train (Hz)
        STEPS = 100   # number of time steps in data
        STEPS_LONG = 5000   # number of time steps in data
        TAU = 0.6      # time constant of calcium indicator (seconds)
        DELTAT = 1/30  # time step duration (seconds)
        self.sigma = 0.1    # standard deviation of gaussian noise
        SEED = 2222    # random number generator seed

        # Make a poisson spike trains
        self.spikes = sima.spikes.get_poisson_spikes(
            deltat=DELTAT, rate=RATE, steps=STEPS, seed=SEED)

        # longer time-series for parameter estimation
        self.spikes_long = sima.spikes.get_poisson_spikes(
            deltat=DELTAT, rate=RATE, steps=STEPS_LONG, seed=SEED)

        # Convolve with kernel to make calcium signal
        np.random.seed(SEED)
        self.gamma = 1 - (DELTAT / TAU)
        CALCIUM = signal.lfilter([1], [1, -self.gamma], self.spikes)
        CALCIUM_LONG = signal.lfilter([1], [1, -self.gamma], self.spikes_long)

        # Make fluorescence traces with random gaussian noise and baseline
        self.fluors = CALCIUM + norm.rvs(
            scale=self.sigma, size=STEPS) + uniform.rvs()
        self.fluors_long = CALCIUM_LONG + norm.rvs(
            scale=self.sigma, size=STEPS_LONG) + uniform.rvs() 

def test_prior():
    """Check that priors are defined for sampling."""
    # define negative log posterior
    posterior_fun = pypesto.Objective(fun=negative_log_posterior)

    # define negative log prior
    prior_fun = pypesto.Objective(fun=negative_log_prior)

    # define pypesto prior object
    prior_object = pypesto.NegLogPriors(objectives=[prior_fun])

    # define pypesto problem using prior object
    test_problem = pypesto.Problem(objective=posterior_fun,
                                   x_priors_defs=prior_object,
                                   lb=-10, ub=10,
                                   x_names=['x'])

    sampler = sample.AdaptiveMetropolisSampler()

    result = sample.sample(test_problem, n_samples=1e4, sampler=sampler,
                           x0=np.array([0.]))

    # get log prior values of first chain
    logprior_trace = -result.sample_result.trace_neglogprior[0, :]

    # check that not all entries are zero
    assert (logprior_trace != 0.).any()

    # get samples of first chain
    samples = result.sample_result.trace_x[0, :, 0]

    # generate ground-truth samples
    rvs = norm.rvs(size=5000, loc=-1., scale=np.sqrt(0.7))

    # check sample distribution agreement with the ground-truth
    statistic, pval = ks_2samp(rvs, samples)
    print(statistic, pval)

    assert statistic < 0.1 

def test_dic_0():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p)
    dics = c.comparison.dic()
    assert len(dics) == 1
    assert dics[0] == 0 

def test_dic_fail_no_posterior():
    d = norm.rvs(size=1000)
    c = ChainConsumer()
    c.add_chain(d, num_eff_data_points=1000, num_free_params=1)
    dics = c.comparison.dic()
    assert len(dics) == 1
    assert dics[0] is None 

def test_bic_data_dependence2():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=2, num_eff_data_points=1000)
    c.add_chain(d, posterior=p, num_free_params=3, num_eff_data_points=500)
    bics = c.comparison.bic()
    assert len(bics) == 2
    assert bics[0] == 0
    expected = 3 * np.log(500) - 2 * np.log(1000)
    assert np.isclose(bics[1], expected, atol=1e-3) 

def test_bic_data_dependence():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=500)
    bics = c.comparison.bic()
    assert len(bics) == 2
    assert bics[1] == 0
    expected = np.log(1000) - np.log(500)
    assert np.isclose(bics[0], expected, atol=1e-3) 

def test_bic_parameter_dependence():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    c.add_chain(d, posterior=p, num_free_params=2, num_eff_data_points=1000)
    bics = c.comparison.bic()
    assert len(bics) == 2
    assert bics[0] == 0
    expected = np.log(1000)
    assert np.isclose(bics[1], expected, atol=1e-3) 

def test_bic_0():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    bics = c.comparison.bic()
    assert len(bics) == 1
    assert bics[0] == 0 

def test_bic_fail_no_num_params():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_eff_data_points=1000)
    bics = c.comparison.bic()
    assert len(bics) == 1
    assert bics[0] is None 

def test_bic_fail_no_data_points():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1)
    bics = c.comparison.bic()
    assert len(bics) == 1
    assert bics[0] is None 

def test_bic_fail_no_posterior():
    d = norm.rvs(size=1000)
    c = ChainConsumer()
    c.add_chain(d, num_eff_data_points=1000, num_free_params=1)
    bics = c.comparison.bic()
    assert len(bics) == 1
    assert bics[0] is None 

def test_aic_data_dependence():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=500)
    aics = c.comparison.aic()
    assert len(aics) == 2
    assert aics[0] == 0
    expected = (2.0 * 1 * 2 / (500 - 1 - 1)) - (2.0 * 1 * 2 / (1000 - 1 - 1))
    assert np.isclose(aics[1], expected, atol=1e-3) 

def test_aic_0():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1, num_eff_data_points=1000)
    aics = c.comparison.aic()
    assert len(aics) == 1
    assert aics[0] == 0 

def test_aic_fail_no_num_params():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_eff_data_points=1000)
    aics = c.comparison.aic()
    assert len(aics) == 1
    assert aics[0] is None 

def test_aic_fail_no_data_points():
    d = norm.rvs(size=1000)
    p = norm.logpdf(d)
    c = ChainConsumer()
    c.add_chain(d, posterior=p, num_free_params=1)
    aics = c.comparison.aic()
    assert len(aics) == 1
    assert aics[0] is None 

def test_aic_fail_no_posterior():
    d = norm.rvs(size=1000)
    c = ChainConsumer()
    c.add_chain(d, num_eff_data_points=1000, num_free_params=1)
    aics = c.comparison.aic()
    assert len(aics) == 1
    assert aics[0] is None 

def evolve(self):
        X = self.state
        dw = norm.rvs(scale=self.dt, size=self.processes)
        dx = self.theta * (self.mu - X) * self.dt + self.sigma * dw
        self.state = X + dx
        return self.state 

def rnorm(n,mean,sd):
    """
    same functions as rnorm in r
    r: rnorm(n, mean=0, sd=1)
    py: rvs(loc=0, scale=1, size=1, random_state=None)
    """
    return norm.rvs(loc=mean,scale=sd,size=n) 

def ou_process_t_driver_batch_fn(num_points, mu, l, sigma, df, x0, dt=1):
    """
    Generates batch of realisation trajectories of Ornshtein-Uhlenbeck process
    driven by t-distributed innovations.

    Args:
        num_points:     int, trajectory length
        mu:             float or array of shape [batch_dim], mean;
        l:              float or array of shape [batch_dim], lambda, mean reversion rate;
        sigma:          float or array of shape [batch_dim], volatility;
        df:             float or array of shape [batch_dim] > 2.0, standart Student-t degrees of freedom param.;
        x0:             float or array of shape [batch_dim], starting point;
        dt:             int, time increment;

    Returns:
        generated data as np.array of shape [batch_dim, num_points]
    """

    n = num_points + 1
    try:
        batch_dim = x0.shape[0]
        x = np.zeros([n, batch_dim])
        x[0, :] = np.squeeze(x0)
    except (AttributeError, IndexError) as e:
        batch_dim = None
        x = np.zeros([n, 1])
        x[0, :] = x0

    for i in range(1, n):
        driver = np.random.standard_t(df, size=df.size) * ((df - 2) / df) ** .5
        # driver = stats.t.rvs(df, loc, scale, size=batch_dim)
        # x_vol = df / (df - 2)
        # driver = (driver - loc) / scale / x_vol**.5
        x[i, :] = x[i - 1, :] * np.exp(-l * dt) + mu * (1 - np.exp(-l * dt)) + \
            sigma * ((1 - np.exp(-2 * l * dt)) / (2 * l)) ** .5 * driver

    return x[1:, :] 

def applyPulseSeq(config, Meq, M0, w, T1, T2, pos0, v, D):
    '''Simulate magnetization vector during nTR (+nDummies) applications of pulse sequence.
    
    Args:
        config: configuration dictionary.
        Meq:    equilibrium magnetization along z axis.
        M0:     initial state of magnetization vector, numpy array of size 3.
        w:      Larmor frequency :math:`2\\pi\\gamma B_0` [kRad/s].
        T1:     longitudinal relaxation time.
        T2:     transverse relaxation time.
        pos0:   position (x,y,z) of magnetization vector at t=0 [m].
        v:      velocity (x,y,z) of spins [mm/s]
        D:      diffusivity (x,y,z) of spins [:math:`mm^2/s`]
        
    Returns:
        magnetization vector over time, numpy array of size [6, nFrames]. 1:3 are magnetization, 4:6 are position

    '''
    M = np.zeros([len(config['t']), 3])
    M[0] = M0 # Initial state

    pos = np.tile(pos0, [len(config['t']), 1]) # initial position
    if np.linalg.norm(D) > 0: # diffusion contribution
        for frame in range(1,len(config['t'])):
            dt = config['t'][frame] - config['t'][frame-1]
            for dim in range(3):
                pos[frame][dim] = pos[frame-1][dim] + norm.rvs(scale=np.sqrt(D[dim]*dt*1e-9)) # standard deviation in meters
            if config['t'][frame]==0: # reset position for t=0
                pos[:frame+1] += np.tile(pos0 - pos[frame], [frame+1, 1])
    if np.linalg.norm(v) > 0: # velocity contribution
        pos += np.outer(config['t'], v) * 1e-6

    for rep in range(-config['nDummies'], config['nTR']): # dummy TRs get negative frame numbers
        TRstartFrame = rep * config['nFramesPerTR']

        for i, event in enumerate(config['events']):
            firstFrame, lastFrame = getEventFrames(config, i)
            firstFrame += TRstartFrame
            lastFrame += TRstartFrame

            M0 = M[firstFrame]

            if 'spoil' in event and event['spoil']: # Spoiler event
                M0 = spoil(M0)

            # frequency due to w plus any gradients 
            # (use position at firstFrame, i.e. approximate no motion during frame)
            wg = w  
            wg += 2*np.pi*gyro*event['Gx']*pos[firstFrame, 0]/1000 # [kRad/s]
            wg += 2*np.pi*gyro*event['Gy']*pos[firstFrame, 1]/1000 # [kRad/s]
            wg += 2*np.pi*gyro*event['Gz']*pos[firstFrame, 2]/1000 # [kRad/s]

            w1 = event['w1'] * np.exp(1j * np.radians(event['phase']))

            t = config['t'][firstFrame:lastFrame+1]
            if len(t)==0:
                raise Exception("Corrupt config['events']")
            M[firstFrame:lastFrame+1] = integrate.odeint(derivs, M0, t, args=(Meq, wg, w1, T1, T2)) # Solve Bloch equation

    return np.concatenate((M, pos),1).transpose() 

def plot_correlation_demo():
    np.random.seed(0)  # to reproduce the data later on
    pylab.clf()
    pylab.figure(num=None, figsize=(8, 8))

    x = np.arange(0, 10, 0.2)

    pylab.subplot(221)
    y = 0.5 * x + norm.rvs(1, scale=.01, size=len(x))
    _plot_correlation_func(x, y)

    pylab.subplot(222)
    y = 0.5 * x + norm.rvs(1, scale=.1, size=len(x))
    _plot_correlation_func(x, y)

    pylab.subplot(223)
    y = 0.5 * x + norm.rvs(1, scale=1, size=len(x))
    _plot_correlation_func(x, y)

    pylab.subplot(224)
    y = norm.rvs(1, scale=10, size=len(x))
    _plot_correlation_func(x, y)

    pylab.autoscale(tight=True)
    pylab.grid(True)

    filename = "corr_demo_1.png"
    pylab.savefig(os.path.join(CHART_DIR, filename), bbox_inches="tight")

    pylab.clf()
    pylab.figure(num=None, figsize=(8, 8))

    x = np.arange(-5, 5, 0.2)

    pylab.subplot(221)
    y = 0.5 * x ** 2 + norm.rvs(1, scale=.01, size=len(x))
    _plot_correlation_func(x, y)

    pylab.subplot(222)
    y = 0.5 * x ** 2 + norm.rvs(1, scale=.1, size=len(x))
    _plot_correlation_func(x, y)

    pylab.subplot(223)
    y = 0.5 * x ** 2 + norm.rvs(1, scale=1, size=len(x))
    _plot_correlation_func(x, y)

    pylab.subplot(224)
    y = 0.5 * x ** 2 + norm.rvs(1, scale=10, size=len(x))
    _plot_correlation_func(x, y)

    pylab.autoscale(tight=True)
    pylab.grid(True)

    filename = "corr_demo_2.png"
    pylab.savefig(os.path.join(CHART_DIR, filename), bbox_inches="tight") 

def plot_mi_demo():
    np.random.seed(0)  # to reproduce the data later on
    pylab.clf()
    pylab.figure(num=None, figsize=(8, 8))

    x = np.arange(0, 10, 0.2)

    pylab.subplot(221)
    y = 0.5 * x + norm.rvs(1, scale=.01, size=len(x))
    _plot_mi_func(x, y)

    pylab.subplot(222)
    y = 0.5 * x + norm.rvs(1, scale=.1, size=len(x))
    _plot_mi_func(x, y)

    pylab.subplot(223)
    y = 0.5 * x + norm.rvs(1, scale=1, size=len(x))
    _plot_mi_func(x, y)

    pylab.subplot(224)
    y = norm.rvs(1, scale=10, size=len(x))
    _plot_mi_func(x, y)

    pylab.autoscale(tight=True)
    pylab.grid(True)

    filename = "mi_demo_1.png"
    pylab.savefig(os.path.join(CHART_DIR, filename), bbox_inches="tight")

    pylab.clf()
    pylab.figure(num=None, figsize=(8, 8))

    x = np.arange(-5, 5, 0.2)

    pylab.subplot(221)
    y = 0.5 * x ** 2 + norm.rvs(1, scale=.01, size=len(x))
    _plot_mi_func(x, y)

    pylab.subplot(222)
    y = 0.5 * x ** 2 + norm.rvs(1, scale=.1, size=len(x))
    _plot_mi_func(x, y)

    pylab.subplot(223)
    y = 0.5 * x ** 2 + norm.rvs(1, scale=1, size=len(x))
    _plot_mi_func(x, y)

    pylab.subplot(224)
    y = 0.5 * x ** 2 + norm.rvs(1, scale=10, size=len(x))
    _plot_mi_func(x, y)

    pylab.autoscale(tight=True)
    pylab.grid(True)

    filename = "mi_demo_2.png"
    pylab.savefig(os.path.join(CHART_DIR, filename), bbox_inches="tight")