<jupyter_start><jupyter_text>Meteoroids planet fall - data synthesisData synthesis as the assignment project for Programming for Data Analysis, GMIT 2019Lecturer: dr <NAME>>Author: **<NAME>**  >Github: [andkoc001](https://github.com/andkoc001/)  >Email: <EMAIL>, <EMAIL>Created: 18-11-2019This Notebook should be read in conjunction with the corresponding README.md file at the project [repository](https://github.com/andkoc001/data_synthesis.git) at GitHub. ___ Introduction  Problem statement and project objectivesPurpose: **To synthesise a dataset that models the phenomenon of meteoroids fall on a planet.**  The phenomenonPhenomenon chosen to simulate in this project is **meteorite planet fall** ([Meteorite](https://en.wikipedia.org/wiki/Meteorite) - Wikipedia, 2019). It is a commonplace event in nature, particularly intensive during planetary systems formation. It occurs when certain category of astronomical objects like asteroids traveling along their trajectory paths hit a planet (or similar celestial objects). As an effect of such an event, an impact with the planet surface often occurs. The crash releases energy, causing resulting in creation of an impact crater.![A new impact crater on Mars](https://www.jpl.nasa.gov/spaceimages/images/wallpaper/PIA22462-800x600.jpg) Image source: [NASA](https://www.jpl.nasa.gov/spaceimages/details.php?id=PIA22462) The scenario and data discussed in the Notebook are fictional, invented solely for the purpose of this project. The data generated here is, however, based on the current understanding of the phenomenon and is deemed to reflect its key characteristics accurately. It could as well simulate a real and true events (for instance in an undiscovered solar system in a galaxy far, far away).Meteoroids planet fall is a highly complex phenomenon than captured in this model. There are a great many factors affecting the phenomenon and relevant to the simulation.  Many of them affect each other, so that the circumstances may appear random and chaotic. Not all the factors were considered in this simulation for simplicity and clarity only. For example: the age and structure of the solar system, the size and density of the planet, the atmosphere thickness, chemical composition of the meteoroids, etc. I believe that the variables chosen for the model reflect the idea of the phenomenon and data synthesised by the simulation reflects its feature adequately.  Glossary of key terms_Complexity_ - ”In science, a complex system is one that is chaotic, and in which the way the system develops feeds back on itself to change the way it is developing.” (Gribbin, 2005)  _Chaos_ - “In science, chaos occurs when a small change in the starting conditions of process produces a big change in the outcome of the process.” (Gribbin, 2005)  _Meteor_ - meteoroids passing the atmosphere and burning in the process.  _Meteorite_ - meteoroid that survived the transit through the atmosphere and hit the ground.  _Meteoroid_ - object of certain size traveling in outer space (beyond the atmosphere of a planet).  _Model_ - a simplified and idealised representation (abstraction) of a system.  _Simulation_ - an approximate imitation of the phenomenon, process or system using its model.    Approach / Method The aim of the project is to create a model of the phenomenon, so that it closely imitates the real event.The first step is identification of the relevant attributes of the phenomenon, based on the subject research. Next, the attributes will be described with variables. Each attribute will require an individual approach. The variables will be modeled by development of algorithms that simulate the real phenomenon characteristics. For this purpose Python capabilities - along with its additional libraries for numerical compuation (NumPy, Scipy, etc) will be employed. Specifically, the NumPy's pseudo-random number generator (PRNG) will be utilised. The results generated through the simulation, will be simultaneously aggregated together into a single dataset.  Simplifications and assumptionsNote: this is a _simplified model_ in terms of factors considered, relationship between them and terminology describing the phenomenon.While designing and developing this simulation, I have made a number of assumptions. Some of them include, inter alia:- A planetary system is subject to celestial mechanics, and is similar enough to our Solar System.- The planet in question is of a rocky type (Earth like) and has an atmosphere (which properties are negligible).- The number of meteorites in the local space is very large, that is running out of the meteorites is not considered.- Units of measurement are arbitrary.   Data categories and typesI have chosen to consider and include in my simulation the following several attributes to describe the phenomenon:1. Time series of observation periods, called in this Notebook as *local* year; in other words - observation sequence number.2. Number of the atmosphere entry, that is an event when a meteoroid gets in contact with the atmosphere  ([Atmospheric entry](https://en.wikipedia.org/wiki/Atmospheric_entry) - Wikipedia, 2019). The meteoroids can then either a) hit the planet, b) burn out completely (due to fraction) while passing the atmosphere, c) skim the atmosphere without hitting the planet surface and continue the fly thorough the space. The two latter were not considered in my simulation.3. Number of the meteorites that survived the transit through the atmosphere (did not burned out completely) and hit the planet surface ([Impact event](https://en.wikipedia.org/wiki/Impact_event) - Wikipedia, 2019).4. Size of the meteorites. This property was idealised in the model and described as mass. Furthermore, only those meteoroids reaching the surface were considered. 5. Relative speed of the meteorites when it hit the planet's surface. Only those hitting the surface were considered.6. Size of the impact crater caused by the falling meteorite - it is a measure of the energy released during the collision; in this simulation I assumed it is a function of only two meteorites parameters: *mass* and *velocity*. NB. only one, the *highest energy*, meteorite in a given period (year) is taken into account for evaluating this variable. *Energy* is calculated using an invented formula that depends on the meteorite's mass and velocity - not recorded in the dataset.7. Classification - meteoroids in my simulation can be of one of the three classes: *iron*, *stony*, *stony-iron*  ([Types of meteorites](https://geology.com/meteorites/meteorite-types-and-classification.shtml) - Geology.com, 2019). NB. only one, the *highest energy*, meteorite in a given year is taken into account for evaluating this variable.These variables will be discussed in greater detail below. ___ Data synthesis  Python environment setup   Importing Python librariesThe following Python libraries were used in synthesising the dataset variables.<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt 
import seaborn as sns
import scipy.stats as stats
# below command will allow for the plots being displayed inside the Notebook, rather than in a separate screen.
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Random number generator seedI have chosen to define seeds for NumPy random methods. That means each time the Notebook and particular cells are run, results will be the same.In order to generate new random results, reset the kernel and run the Notebook cells except for the one below (which sets the seed values).<jupyter_code>determinant = 2019  # seed value pre-set, selected arbitrarily
np.random.seed(determinant)<jupyter_output><empty_output><jupyter_text>Dataset attributesThe above listed attributes were assigned to the following variables in the dataset.* `Year` - observation number.* `Entry` and `Impact` - number of the atmosphere entry and number of hits onto the planet, respectively.* `Mass` and `Velocity` - properties of the meteoroids reaching the planet's surface.* `Crater` - size of the impact crater caused by the falling meteorite* `Class` - classification of the meteoroidNB. units of expression were purposefully not defined in the dataset.Herewith I am creating an empty dataset with headings (attributes) only. The dataset is assigned to the variable `data`.<jupyter_code># creation of empty data (just headings)
data = pd.DataFrame(columns = ["Year", "Entry", "Impact", "Mass", "Velocity", "Crater", "Class"]) 

# data.head() # for testing<jupyter_output><empty_output><jupyter_text>1. Observation number - `Year`The observations are recorded in sequential interval periods of equal length (time series). These periods are time intervals of constant length. Even though I purposefully have not assigned any reference to unit of measurement, it can be interpreted as local years (for the planet in question). This will be assigned to variable **`Year`** in the dataset.Each year is a new observation and a new row in the data table, and is generated by iteration. The total number of observation periods is defined below by an auxiliary variable named `epoch`, therefore the variable defines the size of the dataset.<jupyter_code># defining the length of the epoch, i.e. number of observations
epoch = 1000<jupyter_output><empty_output><jupyter_text>Observation periods will be generated sequentially, starting from year=1, until the end of the epoch. One year means an observation period. The variables will be then assigned to attribute `Year` and incorporated into the dataset.<jupyter_code># assigning year values to the data array, other columns with empty values (NaN - not a number)
year = np.arange(1, epoch+1, 1) # returns evenly spaced values, incrementally modifed by the step value; (start, stop [, step])

# adding the values to the dataset
data["Year"] = year

# data.head(2) # shows first two rows of the current state of the data table; disabled for readibility<jupyter_output><empty_output><jupyter_text>2. Atmosphere entry - `Entry`Number of  meteoroids' entries into the atmosphere of the planet, represented in the dataset by variable **`Entry`**. The intention for this simulation is to generate outputs with relatively low mean and relatively high standard deviation. In reality, the phenomenon is very erratic and chaotic by nature, and can be perceived as a highly unpredictable event. In my model the unpredictability has been reflected by two factors:1. This variable is generated using random method with **[Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution)**. The reason for choosing this distribution is that it describes probability of number of independent events occurring in a time interval. It captures the fact of lesser probability of a larger number of meteoroids entering the planet's atmosphere. The _lambda_ ($\lambda$) value (Poisson distribution parameter) is also generated randomly, to strengthen the uncertainty. The _lambda_ is assigned to an auxiliary variable `entry_lambda`. It will take a value in a half-open range [2, 4).2. An additional amplified random noise. It is assumed that in a given year there may be _none_ to _a very large number_ of entries, with a certain mean number. The noise modifier is assigned to an auxiliary variable `entry_noise`. it will take a value in a range [-1500, -500), chosen by experiments.In order to ensure sensible returns, all negative outputs were converted to positive, using abs() function.The combination of the two sources of randomness, results in the distribution of the generated values taking the shape similar to one that could be expected in a real scenario.<jupyter_code>np.random.seed(determinant) # executing this line ensures the generated numbers are the same

# defines lambda value for the poisson distribution
entry_lambda = np.random.random() + np.random.randint(2, 4) # expected mean value - it gives a general shape of the distribution
print("Lambda:\t", entry_lambda)

# defines noise used to modify the results
entry_noise = np.random.randint(-1500, -500) # extra random noise - expands the range of possible results and adds randomness and increases the standard deviation
print("Noise:\t", entry_noise)
np.random.seed(determinant) # executing this line ensures the generated numbers are the same

# the algorithm with nested random variables - poisson distribution with uniform lambda - and a random noise
# higher probability of occuring fewer meteors entry in a year; only non-negative values
entry = abs(np.random.poisson(abs(entry_lambda + (entry_noise)), len(year)) + entry_noise) # consulted reference: https://stackoverflow.com/a/35730154

# Adding values to the dataset
data["Entry"] = entry<jupyter_output><empty_output><jupyter_text>Analysis of the atmosphere entries<jupyter_code># description of the generated data
#print("atm_entry_lambda:\t", atm_entry_lambda)
#print("atm_entry_noise:\t", atm_entry_noise)
print("Minimum:\t", entry.min())
print("Mean:\t\t", entry.mean())
print("Maximum:\t", entry.max())
print("Std diviation:\t", entry.std())

### plotting - frequency of atmosphere entries
fig, axes = plt.subplots(ncols=2, figsize=(18, 2))
plt.suptitle("Figure 1. Atmosphere entry.", y=1.15, fontsize=14)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=.2, hspace=None)

# left hand plot - number of entries in each year of observation
sns.lineplot(x=year, y=entry, color='royalblue', ax=axes[0])
axes[0].set(xlabel="Year", ylabel="No. of entries", title="Number of atmosphere entries per year.")

# right hand plot - how many the same number of entries occured
sns.countplot(x=entry, ax=axes[1])
axes[1].set(xlabel="Frequency ", ylabel="No. of entries", title="Frequency distribution.")

plt.show()
# Number of meteroids entering the atmosphere, list in sequence of observations
print("\nOccurrences (what numbers happend to be generated):\n", np.unique(set(entry)))

# Assigning the variables
unique_elements, counts_elements = np.unique(entry, return_counts=True) # https://www.w3resource.com/python-exercises/numpy/python-numpy-exercise-94.php

print("\nFrequency of unique values, sorted (Upper value denotess number of occurrences, lower denotes how many time the number occured in all the observations.):\n", np.asarray((unique_elements, counts_elements)))

# print tabularised results - frequency of yearly entries
hits = pd.DataFrame(np.transpose(np.asarray((unique_elements, counts_elements))))
hits.columns = ['Unique values', 'Counter']
hits<jupyter_output>Occurrences (what numbers happend to be generated):
 [{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 64, 65, 66, 67, 68, 69, 70, 71, 72, 75, 76, 81, 83, 92}]

Frequency of unique values, sorted (Upper value denotess number of occurrences, lower denotes how many time the number occured in all the observations.):
 [[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
  24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
  48 49 50 51 52 53 54 55 56 57 58 59 60 61 64 65 66 67 68 69 70 71 72 75
  76 81 83 92]
 [13 21 25 38 29 44 21 27 29 27 30 26 29 30 21 25 25 25 18 24 27 18 21 20
  17 16 20 22 21 15 19 16 18 14 18 12  9 18  8 19 13  9  5  4  9  8  6  3
   3  8  3  4  2  2  4  2  1  4  4  3  5  4  3  2  2  1  1  1  2  1  1  1
   1  1  1  1]]<jupyter_text>3. Surface impact - `Impact`Number of meteorites hits into the planet surface, that is the meteoroids that entered the atmosphere of the planet and penetrated it without being completely burned out in the transit. The variable is represented in my model by variable **`Impact`** which is a result of **binomial distribution** with certain probability of success. The event is simulated by randomly modifying the probability of success, using an auxiliary variable `passing_prob`, which allowable values are in the range [0.05, 0.2) - assumed arbitrarily by experiments. The probability here is generated only once for the entire dataset.This is simplification, as many potentially influencing factors are neglected here that would have an effect in real situation (for instance the meteoroids' size, density, composition, velocity, angle of entry, atmosphere properties).This part is based on dr Ian Mcloughlin's [Notebook](https://github.com/ianmcloughlin/jupyter-teaching-notebooks/blob/master/coin-flip.ipynb) (Mcloughlin, 2019) on coin flip.<jupyter_code>np.random.seed(determinant) # executing this line ensures the generated numbers are the same

# assign probability of meteoroids passing through the atmosfere and hitting the surface
passing_prob = 0.0001*np.random.randint(500, 2001) # range from 0.05 to 0.2

# printing the result
print("Probability of passing through the atmosphere:",round(passing_prob, 4))
np.random.seed(determinant) # executing this line ensures the generated numbers are the same

# filtering the meteoroids - testing how many in each year have passed through the atmosphere

print("Results for a few first observations:\n")

# for each year
for index, row in data.iterrows(): # adapted from https://stackoverflow.com/a/16476974
    #print("\n")
    
    sum = 0 # counter of meteoites successfully passing the atmosphere in a given year
    
    for i in range(data.iat[index, 1]): # data.iat[-1, 1] returns value of a cell
       
        # testing for a single instance (Number of tests, probability of a penetration, number of times to do this);
        passed = np.random.binomial(1, passing_prob, 1) # adapted from https://github.com/ianmcloughlin/jupyter-teaching-notebooks/blob/master/coin-flip.ipynb
        sum = int(sum + passed)
    
    if index <= 5: # print out for first couple of obsevations
        print("Year:", index+1) #, "\tNumber of meteo:\t", row["Entry"])
        print("Number of meteoroids hitting the planet surface / total number entering the atmosphere: ", sum, "/", row["Entry"], "\n")
    
    if row["Entry"] != 0: # check if there was any meteo hittint the atmosphere in a given year
        ratio_passing = sum/row["Entry"]
        
    else:
        ratio_passing = 0  
    
    # assign number of surface impacts in a year to the main dataset
    data.iat[index, 2] = sum<jupyter_output>Results for a few first observations:

Year: 1
Number of meteoroids hitting the planet surface / total number entering the atmosphere:  12 / 37 

Year: 2
Number of meteoroids hitting the planet surface / total number entering the atmosphere:  0 / 5 

Year: 3
Number of meteoroids hitting the planet surface / total number entering the atmosphere:  4 / 32 

Year: 4
Number of meteoroids hitting the planet surface / total number entering the atmosphere:  2 / 12 

Year: 5
Number of meteoroids hitting the planet surface / total number entering the atmosphere:  4 / 32 

Year: 6
Number of meteoroids hitting the planet surface / total number entering the atmosphere:  1 / 8<jupyter_text>The below plot shows number of atmosphere entry against number of successful transits through the atmosphere (and hitting the planet surface).<jupyter_code>data
### plotting - number of impacts to atmosphere entry
fig, axes = plt.subplots(ncols=2, figsize=(18, 2))
plt.suptitle("Figure 2. Surface impacts.", y=1.15, fontsize=14)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=.2, hspace=None)

# left hand plot - number of planet impacts in each year of observation
sns.scatterplot(x="Year", y="Impact", color="darkorchid", alpha=.4, data=data, ax=axes[0])
axes[0].set(xlabel="Year", ylabel="No. of impacts.", title="Number of impacts per year.")

# right hand plot - relationship between no. of atmosphere entries and no. of surface impacts
sns.scatterplot(x="Entry", y="Impact", color="darkorchid", alpha=0.3, data=data, ax=axes[1])   
axes[1].set(xlabel="No. of entries", ylabel="No. of impacts.", title="Relationship between atmosphere entries and surface impacts.")

print("Minimum:\t", data["Impact"].min())
print("Mean:\t\t", data["Impact"].mean())
print("Maximum:\t", data["Impact"].max())
print("Std diviation:\t", data["Impact"].std())
plt.show()<jupyter_output>Minimum:	 0
Mean:		 3.379
Maximum:	 17
Std diviation:	 3.0140410487323233<jupyter_text>4. Mass of the meteorites - `Mass`The mass of the meteorites follows **[Pareto distribution](https://en.wikipedia.org/wiki/Pareto_distribution)**. According to pareto principle, there is many small meteorites and few large. In the dataset this will be assigned to variable **`Mass`**. In the model, the mass of **every** meteorite successfully passing through the atmosphere and hitting the planet surface in given year is generated, but those not reaching the planet's surface are not considered in the dataset. The variable is stored in a data subset for each year (3rd dimension of the original dataset). The variable is generated with the help of `np.random.pareto()` function, which requires passing two parameters:- shape of the distribution - describes how much the number of small meteors would dominate over the large ones; greater the shape value, the more small samples,- mode of the distribution - tells the magnitude (or scale) of the outcome (the size of the meteors); the larger the mode, the larger the scale of the meteors.Below, the pareto function parameters are assigned to auxiliary variables.<jupyter_code># adapted from https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.random.pareto.html

# parameters for pareto dirstribution: pareto_a - shape of the distribution, pareto_m - mode of the distribution

mass_pareto_a, mass_pareto_m = 1.0, 1000.0  # shape and mode - picked arbitrarily to resemble reality as much as the author could think
np.random.seed(determinant) # executing this line ensures the generated numbers are the same

print("Results for a few first observations:\n")

# splicing - getting the number of meteoroids entering the atmosphere in a given year
for i in data["Year"]:
    
    number_meteors = data.loc[i-1]["Impact"] # getting the number of the meteoroids hitting the surface in the subsequnt year
    
    if i <= 5: 
        print("In year", i, "the number of meteoroids passing the atmosphere:", number_meteors)
   
    if number_meteors != 0:  # only those that hit the planet are considered
        
        # genertating mass of each meteoroid
        mass_meteo = np.round((0.01*(np.random.pareto(mass_pareto_a, number_meteors) + 1) * mass_pareto_m) + 0.01*np.random.randint(-1000, 10), 1) # random mass for each meteoroid hitting the surface
        
        # recording the masses in the dataset table
        data.at[i-1, "Mass"] = mass_meteo
        
        if i <= 5:
            print("Their mass are as follows:", data.at[i-1, "Mass"], "\n")
            #sns.distplot(x=number_meteors, y=mass_meteo)#, data=data)
            #sns.scatterplot(x=number_meteors, y=data.at[i-1, "Mass"])
            #print()
                
    else:
        pass<jupyter_output>Results for a few first observations:

In year 1 the number of meteoroids passing the atmosphere: 12
Their mass are as follows: [101.5  14.4  24.5  25.5  81.6  12.2  31.5 101.2  82.2  14.7  16.2  11.6] 

In year 2 the number of meteoroids passing the atmosphere: 0
In year 3 the number of meteoroids passing the atmosphere: 4
Their mass are as follows: [72.3 11.4 32.5 76.9] 

In year 4 the number of meteoroids passing the atmosphere: 2
Their mass are as follows: [20.4 12.3] 

In year 5 the number of meteoroids passing the atmosphere: 4
Their mass are as follows: [10.1 16.2  1.3 26.4]<jupyter_text>The statistical analysis of the meteorites' mass:<jupyter_code>print("Minimum:\t", mass_meteo.min())
print("Mean:\t\t", mass_meteo.mean())
print("Maximum:\t", mass_meteo.max())
print("Std diviation:\t", mass_meteo.std())<jupyter_output>Minimum:	 10.0
Mean:		 60.416666666666664
Maximum:	 294.1
Std diviation:	 104.53193159136697<jupyter_text>In the next cell, there is an analysis of the mass distribution for a larger sample (n=1000) using the above algorithm. The values generated here are separated from the main dataset.<jupyter_code># simulation of the distribution algorithm
mass_distr = np.round((0.01*(np.random.pareto(mass_pareto_a, 1000) + 1) * mass_pareto_m) + 0.01*np.random.randint(-1000, 10), 0)

np.random.seed(determinant) # executing this line ensures the generated numbers are the same

print("Minimum:\t", mass_distr.min())
print("Mean:\t\t", mass_distr.mean())
print("Maximum:\t", mass_distr.max())
print("Std diviation:\t", mass_distr.std())

### plotting - mass distribution
fig, axes = plt.subplots(ncols=2, figsize=(18, 2))
plt.suptitle("Figure 3. Mass distribution pattern (not part of the actual dataset).", y=1.15, fontsize=14)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=.2, hspace=None)

# left hand plot - number of entries in each year of observation
sns.scatterplot(x=year, y=mass_distr, color='darkgreen', alpha=.3, ax=axes[0])
axes[0].set(xlabel="Year", ylabel="Mass", title="Mass of meteroids")

# right hand plot - how many the same number of entries occured
sns.countplot(x=mass_distr, ax=axes[1])   
axes[1].set(xlabel="Mass", ylabel="Frequency", title="Frequency distribution.")

plt.show()<jupyter_output>Minimum:	 7.0
Mean:		 86.955
Maximum:	 12794.0
Std diviation:	 545.4380157038927<jupyter_text>5. Velocity of the meteorites - `Velocity`Meteoroids **`velocity`** variable distribution is an interesting one. Because there is an energy requirement to propel an object, there will be naturally a larger number of slower moving objects, and fewer faster ones. Furthermore, the possible velocity of astronomical objects in space can vary in range from $0$ to $c$ (speed of light). However, only these velocity values can be taken into account which are within the range limited from the bottom by the _minimum velocity_ (planet escape velocity), and from the top - by the the _maximum velocity_. > “On Earth, ignoring the slowing effects of travel through the atmosphere, the lowest impact velocity with an object from space is equal to the gravitational escape velocity of about 11 km/s. The fastest impacts occur at about 72 km/s in the "worst case" scenario in which an object in a retrograde near-parabolic orbit hits Earth. The median impact velocity on Earth is about 20 km/s”> > Source: https://en.wikipedia.org/wiki/Impact_crater In my simulation, for the sake of argument I have chosen the same restriction for the lower and upper limits as for Earth (without units of expression):- _minimum velocity_ = 11,- _maximum velocity_ = 72.For simulating the meteoroids' velocity in the space, I have chosen probability distribution described by a **[normal distribution](https://en.wikipedia.org/wiki/Normal_distribution)** model. The shape of the probability density function was adjusted by manipulating the _expected value_ $\mu$ and _standard deviation_ $\sigma$ to capture the characteristics of the real phenomenon (through experimenting). To capture the fact the phenomenon takes the values only from a restricted domain were simulated (a slice of the probability distribution function ([Truncated distribution](https://en.wikipedia.org/wiki/Truncated_distribution) - Wikipedia, 2019)), I applied the `stats.truncnorm()` function from the SciPy library, which offers to generate the values in a distribution that is truncated from both bottom and top.Similarly to mass variables, values of the velocity are generated for **every** meteorite hitting the planet surface in given observation period - those not reaching the planet's surface are not considered in the dataset. The variable is stored in a data subset for each year (3rd dimension of the original dataset).<jupyter_code>min_velo, max_velo = 11.0, 72.0 # boundary limits - picked arbitrarily to resemble real distribution (as fer as possible)

# parameters for pareto dirstribution: pareto_a - shape of the distribution, pareto_m - mode of the distribution

velo_pareto_a, velo_pareto_m = .3, min_velo  # shape and mode - picked arbitrarily to resemble reality as much as the author could think

# shape of the pareto distribution determines the how much the number of slow meteoroids would dominate over the fast ones; greater the shape value, the more slow samples
# mode in the pareto distribution tells the magnitude (or scale) of the outcome (the velocity of the metoroids); the larger the mode, the larger the velocity of the meteors

### normal distribution properties for stats.truncnorm.rvs() funtion
mu, sigma = -40, 40 # normal distribution parameters
np.set_printoptions(suppress=True, formatter={'float_kind':'{:0.1f}'.format}) # this line eliminates the scientific notation of numbers

np.random.seed(determinant) # executing this line ensures the generated numbers are the same

print("Results for first couple of obsevations:")

# splicing - getting the number of meteors entering the atmosphere in a given year
for i in data["Year"]:
    
    number_meteors = data.loc[i-1]["Impact"] # getting the number of the meteors entereing the atmosphere in the subsequnt year
    
    if i <= 5: 
        print("In year", i, "the number of meteorites passing the atmosphere was:\t", number_meteors) # print out for first couple of obsevations
       
    if number_meteors != 0:  # only those that hit the planet are considered
        
        # generating velocity value for each meteorites
        a = (min_velo - mu) / sigma
        b = (max_velo - mu) / sigma
        velo = np.round(stats.truncnorm.rvs(a, b, loc=mu, scale=sigma, size=number_meteors), 1)
        
        # assigning the velo values to the dataset
        data.at[i-1, "Velocity"] = np.round(velo, 1)
        
        if i <= 5: 
            print("Their velocity was as follows:", data.at[i-1, "Velocity"]) # print out for first couple of obsevations
            print("Min:", velo.min(), "\tMean:", velo.mean(), "\tMax:", velo.max(), "\n")

    else:
         pass<jupyter_output>Results for first couple of obsevations:
In year 1 the number of meteorites passing the atmosphere was:	 12
Their velocity was as follows: [50.2 21.4 30.2 30.8 47.5 18.5 34.1 50.2 47.6 21.8 23.4 17.6]
Min: 17.6 	Mean: 32.775 	Max: 50.2 

In year 2 the number of meteorites passing the atmosphere was:	 0
In year 3 the number of meteorites passing the atmosphere was:	 4
Their velocity was as follows: [14.9 48.5 14.5 65.9]
Min: 14.5 	Mean: 35.95 	Max: 65.9 

In year 4 the number of meteorites passing the atmosphere was:	 2
Their velocity was as follows: [11.7 25.7]
Min: 11.7 	Mean: 18.7 	Max: 25.7 

In year 5 the number of meteorites passing the atmosphere was:	 4
Their velocity was as follows: [15.8 48.1 25.6 28.2]
Min: 15.8 	Mean: 29.425 	Max: 48.1<jupyter_text>Distribution patternHere is a general visual interpretation of the data generated with the above algorithm. To make pattern of the distribution more apparent, number of samples drawn is increased to 1000. We can see, spread of values is significant, but dominate smaller values with only rare occurrences of larger ones.<jupyter_code># simulation of the distribution algorithm
velo_distr = np.round(stats.truncnorm.rvs(a, b, loc=mu, scale=sigma, size=1000), 1)

np.random.seed(determinant) # executing this line ensures the generated numbers are the same

print("Minimum:\t", mass_distr.min())
print("Mean:\t\t", mass_distr.mean())
print("Maximum:\t", mass_distr.max())
print("Std diviation:\t", mass_distr.std())

### plotting - velocity distribution
fig, axes = plt.subplots(ncols=2, figsize=(18, 2))
plt.suptitle("Figure 4. Velocity distribution pattern.", y=1.15, fontsize=14)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=.2, hspace=None)

# left hand plot - number of entries in each year of observation
sns.scatterplot(x=year, y=velo_distr, color='sandybrown', ax=axes[0])
axes[0].set(xlabel="Year", ylabel="Velocity", title="Velocity of meteroids")

# right hand plot - how many the same number of entries occured
sns.violinplot(x=velo_distr, color='sandybrown', ax=axes[1])
axes[1].set(xlabel="Velocity", ylabel="Density", title="Distribution of velocity.")

plt.show()<jupyter_output>Minimum:	 7.0
Mean:		 86.955
Maximum:	 12794.0
Std diviation:	 545.4380157038927<jupyter_text>6. Size of the largest crater - `Crater` In my model, I assume that every meteorite hitting the planet's surface creates an impact crater. However, for simplicity, size of only one crater is simulated in a given year,  caused by the impact of the highest energy meteorite. If there was not impact in a given year, there is also no impact crater. This attribute is designated to variable **`Crater`**.The energy is calculated from the following formula (without specifying the units of measurement):$$ E = \frac {1} {2} m \times v^{2} $$where:$E$ - kinetic energy of the meteorite,$m$ - mass of the meteorite,$v$ - velocity of the meteorite.The calculated value is further randomised by a certain number (noise) that simulate effect of other factors not implicitly considered. The result is always positive by applying the `np.abs()` function.Not every meteorite reaching the planet surface has sufficient energy to create impact crater. This fact in my model is simulated by introducing an threshold value. The crater will be made only if the meteorite's energy is greater than the threshold. The value of threshold is assumed arbitrarily. In the model the impact crater size is directly related (linearly related) to the energy of the meteorite; the size of the crater is factored energy impact (chosen arbitrarily).<jupyter_code># energy threshold value required to create an impact crater
threshold = 80000

# factor between impact energy and crater size caused by the impact
factor = 0.01<jupyter_output><empty_output><jupyter_text>Calculating the impact crater sizes from the highest energy meteorites<jupyter_code>np.set_printoptions(suppress=True, formatter={'float_kind':'{:0.1f}'.format}) # this line eliminates the scientific notation of numbers
np.random.seed(determinant) # executing this line ensures the generated numbers are the same

# effect of random factors - noise
en_noise = np.random.random()

# Calculation of the meteoroids energy
for i in data["Year"]:
    
    if data.loc[i-1]["Impact"] != 0:
        #print(data.loc[i-1]["Mass"], "*", data.loc[i-1]["Velocity"], "=", (data.loc[i-1]["Mass"] * data.loc[i-1]["Velocity"]))
        
        # calculation energy for each meteroid hitting the planet
        energy_calcs = np.round(np.abs(0.5 * data.loc[i-1]["Mass"] * pow(data.loc[i-1]["Velocity"], 2) + en_noise), 1)
        
        # recording only the highest energy metoroid in the dataset (if energy higher than the threshold value)
        if max(energy_calcs) > threshold:
            data.at[i-1, "Crater"] = factor * max(energy_calcs)
        else:
            pass
                        
    else:
         pass<jupyter_output><empty_output><jupyter_text>Statistics of the crater variable.<jupyter_code>print("Frequency:\t", data["Crater"].count())
print("Minimum:\t", data["Crater"].min())
print("Mean:\t\t", data["Crater"].mean())
print("Maximum:\t", data["Crater"].max())
print("Std diviation:\t", data["Crater"].std())<jupyter_output>Frequency:	 170
Minimum:	 812.582
Mean:		 5378.922194117646
Maximum:	 99210.03800000002
Std diviation:	 13155.850198023258<jupyter_text>7. Class of the highest energy meteoroid Even though the relevant science recognises a multi-level meteorites taxonomy ([Types of meteorites](https://geology.com/meteorites/meteorite-types-and-classification.shtml) - Geology.com, 2019; [Meteorite classification] (https://en.wikipedia.org/wiki/Meteorite_classification) - Wikipedia, 2019), I assumed in my model the meteorites can be one of the following three classes only based on their material.* iron,* stony,* stony-iron.Each class has a certain probability of of occurrence - see below.<jupyter_code># in the solar system occure the following classes of meteorids as well as their corresponding probability of occurance
classes, class_probability = ("iron", "stony", "stony-iron"), (0.86, 0.12, 0.02) # tuple data types, so that the order matters
np.random.seed(determinant) # executing this line ensures the generated numbers are the same

for i in data["Year"]:
    
    if data.loc[i-1]["Impact"] != 0: # on condition there was an impact in this year
        
        # assigning the class in given year
        data.at[i-1, "Class"] = np.random.choice(classes, p=class_probability)
        
    else:
         pass<jupyter_output><empty_output><jupyter_text>Here is distribution of the classes that fall caused creation of the impact crater.<jupyter_code>sns.catplot(x="Crater", y="Class", height=2, aspect=6, data=data).fig.suptitle("Figure 5. Craters made by each class of the meteorites.", y=1.0, fontsize=14)
plt.show()<jupyter_output><empty_output><jupyter_text>___ Final datasetThe final form of the synthesised dataset is shown below.<jupyter_code>data<jupyter_output><empty_output><jupyter_text>Visualisation of the entire dataset<jupyter_code># Pair plot of the final dataset
sns.pairplot(data, hue="Class").fig.suptitle("Figure 6. Final dataset.", y=1.0, fontsize=14)
plt.show()<jupyter_output>/home/ak/anaconda3/lib/python3.7/site-packages/statsmodels/nonparametric/kde.py:447: RuntimeWarning: invalid value encountered in greater
  X = X[np.logical_and(X > clip[0], X < clip[1])] # won't work for two columns.
/home/ak/anaconda3/lib/python3.7/site-packages/statsmodels/nonparametric/kde.py:447: RuntimeWarning: invalid value encountered in less
  X = X[np.logical_and(X > clip[0], X < clip[1])] # won't work for two columns.
{'fraction_non_alphanumeric': 0.08541716281654362, 'fraction_numerical': 0.03183297610542193, 'mean_word_length': 4.569933677229182, 'pattern_counts': {'":': 0, '<': 86, '<?xml version=': 0, '>': 88, 'https://': 21, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '624090', 'license': 'MIT', 'n_tokens_mistral': 11664, 'n_tokens_neox': 10281, 'n_words': 4767, 'path': 'data_synth_1.ipynb', 'repo_name': 'andkoc001/data_synthesis', 'size': '37833'}