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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
|
<html>
<head>
<title>Probability Distributions</title>
</head>
<body>
<h1>Probability Distributions</h1>
<p>
Open the file <tt>examples/cain/CaoPetzold2006_Schlogl.xml</tt>.
(Note that this model is used in the
<a href="VisualizationPlottingHistograms.htm">Plotting Histograms</a>
section of the
<a href="Visualization.htm">Visualization and Analysis</a>
chapter.)
For the "Time Series" method, change the number of frames
to 6. Thus, we record the state every two time units.
Generate 100 trajectories and then click the export button
<img src="filesave.png"> in the simulation output panel.
In the export configuration window, deselect the
"Time in first column" checkbox. Then export
the time series data with the base name "schlogl".
</p>
<p>
In R, we first read the data into a data frame, which has 6 rows and
100 columns. We extract the row for time t = 10 as a matrix.
Then we plot a histogram of the data.
</p>
<pre>
> schlogl <- read.csv('schlogl.csv')
> dim(schlogl)
[1] 6 100
> x <- as.matrix(schlogl[6,])
> hist(x)
</pre>
<p align="center">
<!--Save as 8x6. Export to PNG and then resize to 4x3.-->
<img src="ExportR/Probability/hist.png">
</p>
<p>
Next we customize the histogram. We increase the number of breaks,
which increases the number of bins. We switch from frequency to
probability density. Finally, we specify a main title and the y axis
label.
</p>
<pre>
> hist(x, breaks=31, freq=F, main='Population at t = 10', xlab='Population')
</pre>
<p align="center">
<img src="ExportR/Probability/histCust.png">
</p>
<p>
Another possibility for visualizing the probability distribution is to
use kernel density estimation. This gives a smoothed version of
the distribution. We plot the density and then use the <tt>rug()</tt>
function to show the data points below the graph.
</p>
<pre>
> plot(density(x), main='Probabilty density at t = 10', xlab='Population')
> rug(jitter(x))
</pre>
<p align="center">
<img src="ExportR/Probability/density.png">
</p>
<p>
Note that you may directly export histogram data from Cain. However,
when working with R, it is usually better to export the raw time
series data.
</p>
</body>
</html>
|
