<jupyter_start><jupyter_text>Stochastic Processes ExercisesIn this notebook we'll explore stochastic processes. Previously, we looked at a periodic signal, which is purely *deterministic*. This implies that for any given point in time $t_i$, you can exactly compute the value of $y_i$. For stochastic processes, this is not true. For any time $t_i$, you can *at best* compute a probability distribution for $y_i$. There are several degrees of randomness. A pure white noise process has data points that are $strictly independent$. Each data point is sampled from some distribution $p(y_i)$, which does not pretend on any of the previous time points.As explained above, deterministic processes are those where any time point is exactly and unique determined by some functional form, $y_i = f(t_i)$. In between lies the realm of correlated stochastic processes and dependent, but uncorrelated processes. For the former, there is no unique functional form $f(t_i)$ that would tell us what $y_i$ is, but neither is $y_i$ determined by $p(y_i)$. Instead, in order to calculate the current value of $y_i$ might be determined by some probability distribution $p(y_i)$, but could at the same time also be a function of $g(t_{i-1})$, $h(t_{i-2})$ and so on.  A running sum of coin tossesIn the lecture, we've looked at an example that involved a running sum of coin tosses.The rules are as follows:* start at `cointoss_sum=0`* at every step, toss a coin (it can be a virtual one)* if the coin comes up heads, add $1$ to the running sum* if the coin comes up tails, subtract $1$ from the running sum* Do this experiment 10000 times and plot the resulting time series (Hint: `numpy.random.choice` is very helpful here) The coin toss experiment is very close to the *autoregressive (1) process* (AR(1) process). Autoregressive processes are processes where the current data points depends on the $n$ previous data points. For the AR(1) process, the current data point only depends on one previous data point:    $$ y_i = c + \phi y_{i-1} + \epsilon_{\sigma_n} $$    where $c$ is a constant, $\phi is a parameter and $epsilon_{\sigma_n}$ is a normally distributed random variable with a mean of $0$ and a standard deviation of $\sigma_n$. * Define a function that creates an AR(1) process and create a time series with $c=1$, $\phi=0.5$ and $\sigma_n = 0.1$. Plot the result * Is this process stationary? Hint: Make a long time series, subdivide into segments and compute the mean and variance of each. Plot the resulting distributions.* (Bonus) Also look at the autocorrelation function for the various segments. Does the ACF depend on time? * Is the process ergodic? Why (not)?* Could the time series you've generated be a realistic light curve from an astronomical object? Why (not)? Astronomical time series cannot go below zero. Thus, in order to make this process into a time series that could represent brightness, we need to make additional assumptions. In many sources (e.g. AGN and X-ray binaries), a good assumption, backed by observations, is that what we observe is actually the exponential of a stochastic process like the one defined above. * Make a time series from the AR(1) process with 100000 samples. Exponentiate the resulting time series and plot the result. What does this time series look like to you? * Compute the power spectrum (in Leahy normalization) of the original time series, as well as the exponentiated time series. Has exponentiation changed the shape of the power spectrum? (Hint: you're going to have to choose a time scale and make a corresponding list of time stamps): In realistic applications, measurements almost always include some kind of detector noise component. Here, we've already included Gaussian noise (in the definition of the AR(1) process above. You could replace that with Poisson noise as in yesterday's exercises.  Averaging Power SpectraIn the lecture, we've learned that power spectra of single time series are very noisy. One useful way to reduce the noise in the power spectrum is to average lots of individual segments. This works for most stochastic processes.* Make a long observation of your AR(1) process with many data points with a constant rate and Gaussian uncertainties. Exponentiate the data set.* Segment your light curve into intervals of equal duration, compute the power spectrum for each and then average the powers. * Plot and compare the resulting power spectrum with that of the complete, unsegmented light curve* Try out different lengths for the segment intervals, such that you average different numbers of segments (e.g. [2, 5, 10, 20, 50, 100, 500]). What trade-off do you need to make when averaging power spectra? (Hint: look at both the uncertainty distribution of the powers as well as the frequencies)* Stingray has an `AveragedPowerspectrum` class that helps you do much of the work above automatically. The syntax is quite similar to that of `Powerspectrum`. Look at the docstring and make an averaged power spectrum using that class.  Parametric modeling of stochastic processesIn practice, we'd like to model the power spectrum with a parametric model so we can estimate the parameters (e.g. the power law index), which might directly relate to physical processes. For this, we're going to have to write down a likelihood and some priors.For now, let's assume that our model is a  power law and a constant (for the white noise level at high frequencies):* Write down a function for a power law (with a negative exponent) and a constant . The function should take as input the x-variable (frequency) and the three parameters (power law index, amplitude and constant white noise level).* Make a set of frequencies going from $0.1$ Hz to $100$ Hz. Make a power law with a power law index of -2, an amplitude of 10 and a white noise level of 2* plot the result For this exercise, we're going to use the coin toss experiment rather than the AR(1) process, because it has a much simpler power spectrum we can actually approximate by a power law.* Make a time series of the coin toss experiment, where the time array is 16 seconds long and there are 1024 data points per second. Note: instead of using [1, -1], use [0.04, -0.04] to reduce the overall variance in the time series.* exponentiate the resulting time series* for each data point, pick from a Poisson distribution with the rate parameter equal to the value of the exponentiated time series* make a power spectrum out of that data set and plot the result In order to model our (averaged) power spectrum above, we now need to define likelihood, prior and posterior. The posterior should only take a list of parameters as input.As we've learned in the lecture (and the previous exercise), the power spectrum is distributed as $\chi^2$ with $m$ degrees of freedom, where $m$ is the number of segments you've averaged. In practice, the log-likelihood for this distribution is defined as (see [<NAME>, 2012](http://adsabs.harvard.edu/abs/2012ApJ...746..131B) for details): $$\log{ p(\mathrm{Data} | \theta, I)} = -M \sum{\{\frac{P_j}{S_j} + \ln{S_j} + \left( \frac{1}{M} - 1 \right) \ln{P_j} + C(2M) \}}$$where $P_j$ is the power at frequency $f_j$, $S_j = g(f_j, \theta)$ is the value of the model for the power spectrum $g$ with parameters $\theta$ at frequency $f_j$, $M$ is the number of averaged segments (1 if you haven't averaged any segments). $C(2M)$ is a constant for a given $M$, so for parameter estimation, we can simply ignore it. * In many practical applications, astronomers use the standard Gaussian likelihood instead of the $\chi^2$ likelihood above. Do you think this is a good idea? Why (not)? Can you think of circumstances where that might be appropriate?* Make a python function of the $\chi^2$ likelihood for 2 degrees of freedom. Hint: Calculate the natural logarithm instead of the likelihood itself. Your likelihood needs to    * compute the model for the given parameter    * use the model and the data to compute the likelihood    * return the likelihood* Use the original power spectrum (not averaged) and compute the likelihood of that model for a range of parameters (play around with a power law index between [0, 4], an amplitude between [0.001, 1000] and a noise level between [1.0, 3.0])* Which parameters give you a large likelihood?Hint: sometimes, the log-likelihood might become `inf` or `nan` for exceptionally strange values. The standard way to deal with this is to catch this case and return a very small number instead.<jupyter_code># very small value to substitute of prior is inf or nan
logmin = -1.0e16

def loglikelihood(pars):

    # compute model power spectrum
    ps_model = 
    
    # compute log-likelihood
    loglike = 
    
    # catch cases where the log-likelihood is either inf or nan:
    if not np.isfinite(loglike):
        return logmin
    else:
        return loglike<jupyter_output><empty_output><jupyter_text>Exploring parameters by hand is time-consuming and no fun. There are better methods. You can use either *optimization* to find the maximum of the distribution (also called the mode) or *sampling* to characterize the posterior distribution.For this, however, we need to define the posterior distribution instead of just the likelihood.* What do you think would make good priors for this model? Think back to the earlier lectures this week (if you can't figure this out, some suggestions are further down below). Hint: The power law parameter can not be negative, and also not very large (probably  not more than 10 for this model); the amplitude, conversely could span several orders of magnitude; assuming you're computing power spectra in Leahy normalization, your constant noise level should be fairly tightly clustered around a value of $2$.* define a `logprior` function that takes a list of the parameters as input and returns the logarithm of the prior. Hint: for uniform priors with boundaries, the log-prior might become `inf` or `nan`. The standard way to deal with this is to catch this case and return a very small number instead. Here are some suggestions for priors for the three parameters in the model:* power law index: uniform prior between (0,5)* amplitude: a Jeffrey's prior (i.e. a uniform prior on the log(amplitude)) with wide boundaries (say, (-20, 20) in log(amplitude) space)* background noise level: a normal distribution with a mean of $2$ and a fairly narrow standard deviation (maybe $0.1$). Note: a Jeffrey's prior on the log(amplitude) means that your posterior should actually take the log(amplitude) as a parameter rather than the amplitude itself. This means you might need to change your likelihood above to exponentiate the amplitude before computing the model power spectrum.* If you haven't already, now write a function for the log-prior and the log-posterior. Hint: the log-posterior is essentially the sum of the log-prior and the log-likelihood.<jupyter_code>def logprior(pars):
    index = pars[0]
    # uniform prior for the power law index between 0 and 5
    # make index 1/5 inside that range, and zero everywhere else
    p_index =
    
    # exponential prior for the log-amplitude
    logamplitude = pars[1]
    # uniform prior on log-amplitude between -20 and 20
    p_logamplitude = 
    
    whitenoise = pars[2]
    # define a normal distribution with a mean of 2 and a standard deviation of 0.1
    norm_dist = 
    p_whitenoise = 
    
    logprior = 
    
    if not np.isfinite(logprior):
        # add code here
        
    else:
        return logprior<jupyter_output><empty_output><jupyter_text>Okay, cool! We're now ready to work with that. * Compute the posterior for a few different parameters. How does it compare to the likelihood?<jupyter_code>import scipy.optimize
res = scipy.optimize.minimize(mylikelihood, x0, method='L-BFGS-B')<jupyter_output><empty_output><jupyter_text>OptimizationLet's start with optimization. Here, we essentially fit the model to the data such that we find the maximum of the posterior distribution. In general, optimization is a hugely complex topic, and no one algorithm works perfectly for all applications. IN general, however, we're going to get fairly decent results for fitting power spectra with some of the standard algorithms.Optimization algorithms are implemented in `scipy.optimize`. One particular quirk of these algorithms is that they find the *minimum* of a function rather than the maximum. Since we're interested in the maximum of the posterior, we'll need to *invert* the posterior in order to get a result. Do that below. Okay, cool. Now we're ready to use the optimizer. I've written that down for you, since I don't want you to spend your exercise time digging through the `scipy.optimize` [documentation](http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html). It can be really useful, though, in case you'd like to try out a different algorithm or different parameters for the algorithm you've chosen. I've picked the `L-BFGS-B` algorithm because it's generally quite robust.You're going to have to guess some parameters in `x0` to start the optimization.<jupyter_code>import scipy.optimize
res = scipy.optimize.minimize(inveres_posterior, x0, method='L-BFGS-B')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.061286178493161036, 'fraction_numerical': 0.010428474268873271, 'mean_word_length': 4.428219852337982, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '397577', 'license': 'MIT', 'n_tokens_mistral': 3453, 'n_tokens_neox': 3235, 'n_words': 2006, 'path': 'imprs_stochasticprocesses.ipynb', 'repo_name': 'dhuppenkothen/imprs_materials', 'size': '13245'}