These days more and more rockets are developed by companies, students and other amateurs. To assess the safety of these rockets, a method is necessary to estimate the risk of these vehicles quantitatively. Because tools to do so are not freely available, the goal of this thesis is to research methods to quantitatively estimate the risk of nominal flight of unguided rockets. The first step to do so, is to develop a trajectory simulator. To incorporate the effect of wind (weather cocking) and spinning, it is necessary to develop a simulator with three translational degrees of freedom, and three rotational degrees of freedom. The position and velocity of the rocket are described in a Cartesian coordinate system. To prevent singularities in the orientation, quaternions are used to describe the orientation of the rocket. To propagate these through time, classical equations of motion for a rigid body can be used. To do so the impulse thrust force and the jet damping moment are considered as an external force and external moment. The gravity force, the aerodynamic force, and the aerodynamic moment are the other external influences. To calculate these, the Earth is modeled as a flatted sphere (both shape and gravity field) with a US 1976 standard atmosphere, and a tabulated wind profile. The aerodynamic forces and moments are calculated based on this environment and a limited set of aerodynamic coefficients, obtained from the Missile Datcom 97 program. Because the problem’s nature it is not solvable analytically. Furthermore, in some parts of the flight the state rocket’s state changes rapidly, while in others parts the system behaves very smoothly. Therefore a Runge-Kutta-Fehlberg 56 variable time step integrator is used to solve the equations of motion numerically. As a result of this decision it is necessary to implement an event handling algorithm to execute discrete behavior (e.g. stage ignition) at the right time. To validate the in C++ developed program ROSIE, a small amateur rocket’s flight (2m long, apogee at 1 km) is simulated. The results show similar behavior as observed in the real flight data of the simulated rocket, indicating that the simulator is valid for the simulation of small amateur rockets. When comparing a single simulated impact point with a measured impact point, they will most likely not coincide. This can be explained by errors in the flight’s modeling, and can be due to randomness in what is simulated. By stating that these errors are due to uncertainties in the analysis, and modeling these uncertainties as random (input) variables, the result of the analysis will also be a random variable. An advantage of such a variable is that the result is not a single point, but a probability distribution. This distribution can be used to determine the risk to each individual, or the area which needs to be evacuated. To determine the probability distribution different methods are used. These methods can be grouped in parameterized and non-parameterized methods. With parameterized methods, a known type of probability distribution is used to model the dispersion of the impact point. To estimate the parameters of this distribution only a limited number of simulation runs is necessary. The downside of using a parameterized method is that the problem is simplified, in such a way that it can be modeled by a known distribution function. This introduces errors in the estimation of the risk. The alternative is to use a non-parameterized method. The method used is Monte Carlo simulation. To determine the probability of impact in a certain area, the trajectory simulation is run a large number of times. Because part of the input is random, each estimated impact point will have another location. By counting the number of impact points which landed in the sub domain, and dividing this by the total number of simulations, the probability of impact in that specific area can be calculated. Properties of Monte Carlo simulation can be used to specify a confidence interval for the calculated answer. By increasing the number of simulations the accuracy of the answer can be increased. The downside of Monte Carlo simulation is that many trajectory simulations are necessary, making it a very computational time intensive method. For the parameterized methods two different types of probability distributions are used, being the bivariate normal distribution, and the mixture model variant. This variant combines multiple bivariate normal distribution into a new more complex model. To estimate the parameters of these distributions, they are fitted through the data generated with trajectory simulator. This data is generated with different sampling schemes. Four different schemes are used. One-at-a-time sampling and orthogonal array sampling are only valid for the bivariate normal distribution, while the random sampling and Sobol sequence sampling scheme can be used for both. To compare the different methods, they are applied to a preliminary design of the Stratos II rocket. This is a large amateur rocket (5m) which will fly the 50 km in the near future. As a starting position for the comparison, 50000 random trajectories are simulated as benchmark data (Monte Carlo simulation). Comparing the parameters of the bivariate normal distribution obtained with this simulation and with 85 one-at-a-time samples shows a differences of 15 to 30 %. When 81 orthogonal array samples are used, this difference is reduced to less than 10%. With the same number of samples the found accuracy of Sobol sequence sampling is similar to that of the orthogonal arrays, but an advantage of Sobol sampling is that any number of samples can be selected to reach the desired accuracy. Random sampling cannot be used at a low number of samples. Orthogonal array and Sobol sampling have a better performance than one-at-a-time sampling, but are also more complex in their use. Which method to use, depends on the application. The second aspect which is benchmarked, is how accurate the parameterized method estimated the probability of impact at two specific locations. Comparing the probability (and its confidence interval) obtained directly from the Monte Carlo simulation with the one obtained from the bivariate normal distribution (based on 50000 random points) shows that the latter method can be either a factor 2 to 3 too conservative, or too optimistic. Whether or not this is acceptable depends on the applications. Using a mixture model, consisting of only two bivariate normal distributions, already shows significant improvements. Increasing the number of components does increase the accuracy, but also the sensitivity to outliers. Therefore it is not wise to use more than 5 components. To determine the influence of each individual parameter on the probability distribution (the footprint), different sensitivity methods are applied. To study non-linear effects and interactions, the set of parameters is reduced with the Morris elementary effect screening method. Applied regression based sensitivity analysis shows that the dispersion in range direction consists mainly of linear contributions. More specifically, these are the contributions due to uncertainty in thrust of the second stage, the drag coefficient and the atmospheric density. The dispersion in cross range direction depends for a significant part (> 10%) on non-linear contributions and interactions between parameters. This explains the error when modeling the probability distribution with a bivariate normal distribution. For the studied case, the two parameters which have the most influence on the dispersion in cross range direction, are the thrust misalignment of the first stage and the center of mass offset. The Stratos II is simulated to be launched from Esrange Space Center. Comparing the found footprint with the allowed impact area of Esrange shows that it is unlikely that the studied design is allowed to be launched from this site. By improving the design, it might be possible to reduce the size of the footprint in such an extend that it will fit if nominal flight occurs. Unfortunately, due to the experimental nature of the rocket, Esrange will likely demand that for non-nominal flight too, the footprint has to be in the allowed impact area. It is very unlikely that this is possible for a rocket flying to 50 km.