{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pseudo random number generation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate two sequences of 10 uniform deviatives using numpy.random.uniform() and plot them"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"samp1=\n",
"samp2=\n",
"plt.plot(samp1,'rs')\n",
"plt.plot(samp2,'go')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now demonstrate how you can generate a reproducible set of \"random\" numbers by setting a seed using numpy.random.seed() before each of your calls:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# set a seed first\n",
"\n",
"samp1=\n",
"# set the seed again to the same value\n",
"\n",
"samp2=\n",
"plt.plot(samp1,'rs')\n",
"plt.plot(samp2,'go')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"

### Generating random deviates from a non-standard probability distribution function

\n",
"\n",
"Suppose we want to generate random numbers for a \"right triangular\" PDF, i.e. $P(x) \\propto x$ for $x$ between 0 and 1, and zero everywhere else"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot to show the \"right triangular\" PDF\n",
"x=[-1,0,1,1,2]\n",
"y=[0,0,1,0,0]\n",
"plt.plot(x,y)\n",
"ticks=plt.yticks(y,' ')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, what is expression for the properly normalized PDF, i.e. a function that looks like this which integrates to unity?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**ANSWER HERE:** "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will generate the random deviates using the transformation method. What is the expression for the cumulative distribution function $CDF(x)$?\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**ANSWER HERE:** "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What is the expression for the inverse of the CDF, i.e. given $0 ANSWER HERE: "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given this, write a routine to generate a random deviate in the desired function: generate a uniform deviate in the CDF, and transform it back to $x$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def triangle(size=1) :\n",
" \"\"\" Generate size random deviates in a triangular PDF\"\"\"\n",
" \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use your routine to generate samples with different sizes and plot the histograms to validate your routine"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#generate sample using your routine and plot histogram\n",
"samp=triangle(size=100000)\n",
"plt.hist(samp,bins=np.arange(0,1,0.01))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"OK, how about a triangular PDF, i.e one that is proportional to x for 0 ANSWER HERE: \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Implement code to generate random deviates from this PDF."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def triangle(size=1) :\n",
" \"\"\" routine to return random deviates in a triangular function\"\"\"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#generate and plot your deviates\n",
"samp=triangle(size= ) # fill in sample size\n",
"plt.hist(samp,bins ) # fill in bin sizes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"OK, let's do an astronomical example. Stars come in a range of masses, with low mass stars being more common than high mass stars. A simple parameterization of the so-called initial mass function is the Salpeter mass function, which has the form $P(M)\\propto M^{-2.35}$ for $0.1 < M < 100$, where $P(M)$ is the probability of getting a star with mass between M and M+dM. \n",
"\n",
"Write a routine to generate samples from a power law distribution, e.g., for an initial mass function of the form $P(M)\\propto M^{-2.35}$ for $0.1 < M < 100$ (a Salpeter IMF). One might use this, for example, to study how well one could constraint the power law slope of a stellar mass function, given some finite number of stars (e.g., how many stars do you need to measure the IMF slope to an accuracy of XX?).\n",
"\n",
"Start by deriving the expressions for the normalized PDF, the CDF, and the inverse of the CDF"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"** ANSWER HERE: **"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now implement this in some routines"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def pdf(x) :\n",
" \"\"\"routine to return normalized continuous Salpeter IMF given input mass\"\"\" \n",
" \n",
"x=np.arange(0.1,100,0.1)\n",
"plt.loglog(x,p(x))\n",
"def salpeter(size=1) :\n",
" \"\"\" routine to generate a sample drawn from Salpeter IMF of input size\"\"\"\n",
" \n",
" \n",
"samp=salpeter(size= )\n",
"plt.figure()\n",
"plt.hist(samp,bins= )\n",
"plt.xscale('log')\n",
"plt.yscale('log')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Does your sample match the desired mass function?\n",
"

** ANSWER HERE **:\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Additional examples:

\n",
"\n",
"It is generally thought that the Salpeter IMF has too many faint stars. More realistic IMFs include those characterized by Kroupa et al, and by Chabrier et al. See plot .\n",
"\n",
"The Kroupa et al IMF is characterized by a 3-piece power law:

\n",
" $dN/dM\\propto M^{-2.7}$ for $M>1$ $M_{sun}$

\n",
" $dN/dM\\propto M^{-2.2}$ for $0.5\n",
" $dN/dM\\propto M^{-1.3}$ for $M<0.5$ $M_{sun}$

\n",
"\n",
"The Chabrier IMF is charaterized by a lognormal distribution at masses less than a solar mass, and a power law (e.g. with Salpeter slope) above it:

\n",
"$dN/dM \\propto \\exp -(\\log M - \\log M_0)^2$ for $M Schechter luminosity function , which declines exponentially towards high luminosities, and goes as a power law towards low luminosities:

\n",
"$dN/dL \\propto \\exp{-L/L_*} L^\\alpha$\n",
"\n",
"Or you could consider some other distribution related to your research....\n",
"\n",
"If you were so motivated, you could write routines to generate samples from one of these. Or you could describe how you would go about doing so....:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.18"
}
},
"nbformat": 4,
"nbformat_minor": 4
}