<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Salmon
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License
Step1:  Can we predict salmon populations?
Each year the U.S. Atlantic Salmon Assessment Committee reports estimates of salmon populations in oceans and rivers in the northeastern United States.  The reports are useful for monitoring changes in these populations, but they generally do not include predictions.
The goal of this case study is to model year-to-year changes in population, evaluate how predictable these changes are, and estimate the probability that a particular population will increase or decrease in the next 10 years.
As an example, I'll use data from page 18 of the 2017 report, which provides population estimates for the Narraguagus and Sheepscot Rivers in Maine.
There are tools for extracting data from a PDF document automatically, but for this example I will keep it simple and type it in.
Here are the population estimates for the Narraguagus River
Step2:  To get this data into a Pandas Series, I'll also make a range of years to use as an index.
Step3:  And here's the series.
Step4:  Here's what it looks like
Step5:  Modeling changes
To see how the population changes from year-to-year, I'll use diff to compute the absolute difference between each year and the next.
Step6:  We can compute relative differences by dividing by the original series elementwise.
Step7:  These relative differences are observed annual net growth rates.  So let's drop the 0 and save them.
Step8:  A simple way to model this system is to draw a random value from this series of observed rates each year.  We can use the NumPy function choice to make a random choice from a series.
Step9:  Simulation
Now we can simulate the system by drawing random growth rates from the series of observed rates.
I'll start the simulation in 2015.
Step10:  I'll create a System object with variables t_0, p_0, rates, and duration=10 years. 
The series of observed rates is one big parameter of the model.
Step11:  Write an update functon that takes as parameters pop, t, and system.
It should choose a random growth rate, compute the change in population, and return the new population.
Step12:  Test your update function and run it a few times
Step14:  Here's a version of run_simulation that stores the results in a TimeSeries and returns it.
Step15:  Use run_simulation to run generate a prediction for the next 10 years.
The plot your prediction along with the original data.  Your prediction should pick up where the data leave off.
Step17:  To get a sense of how much the results vary, we can run the model several times and plot all of the results.
Step18:  The plot option alpha=0.1 makes the lines semi-transparent, so they are darker where they overlap.
Run plot_many_simulations with your update function and iters=30.  Also plot the original data.
Step20:  The results are highly variable
Step21:  Test your function by running it with iters=5.
Step22:  Now we can run 1000 simulations and describe the distribution of the results.
Step23:  If we substract off the initial population, we get the distribution of changes.
Step24:  The median is negative, which indicates that the population decreases more often than it increases.
We can be more specific by counting the number of runs where net_changes is positive.
Step25:  Or we can use mean to compute the fraction of runs where net_changes is positive.
Step26:  And here's the fraction where it's negative.
Step27:  So, based on observed past changes, this model predicts that the population is more likely to decrease than increase over the next 10 years, by about 2
Step28:  I'll add the weights to the System object, since they are parameters of the model.
Step29:  We can pass these weights as a parameter to np.random.choice (see the documentation)
Step30:  Write an update function that takes the weights into account.
Step31:  Use plot_many_simulations to plot the results.
Step32:  Use run_many_simulations to collect the results and describe to summarize the distribution of net changes.
Step33:  Does the refined model have much effect on the probability of population decline?
<ASSISTANT_TASK:>
Python Code:
# install Pint if necessary
try:
    import pint
except ImportError:
    !pip install pint
# download modsim.py if necessary
from os.path import exists
filename = 'modsim.py'
if not exists(filename):
    from urllib.request import urlretrieve
    url = 'https://raw.githubusercontent.com/AllenDowney/ModSim/main/'
    local, _ = urlretrieve(url+filename, filename)
    print('Downloaded ' + local)
# import functions from modsim
from modsim import *
Explanation: Salmon
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International
End of explanation
pops = [2749, 2845, 4247, 1843, 2562, 1774, 1201, 1284, 1287, 
        2339, 1177, 962, 1176, 2149, 1404, 969, 1237, 1615, 1201]
Explanation: Can we predict salmon populations?
Each year the U.S. Atlantic Salmon Assessment Committee reports estimates of salmon populations in oceans and rivers in the northeastern United States.  The reports are useful for monitoring changes in these populations, but they generally do not include predictions.
The goal of this case study is to model year-to-year changes in population, evaluate how predictable these changes are, and estimate the probability that a particular population will increase or decrease in the next 10 years.
As an example, I'll use data from page 18 of the 2017 report, which provides population estimates for the Narraguagus and Sheepscot Rivers in Maine.
There are tools for extracting data from a PDF document automatically, but for this example I will keep it simple and type it in.
Here are the population estimates for the Narraguagus River:
End of explanation
years = linrange(1997, 2015)
years
Explanation: To get this data into a Pandas Series, I'll also make a range of years to use as an index.
End of explanation
pop_series = TimeSeries(pops, index=years)
pop_series
Explanation: And here's the series.
End of explanation
def plot_population(series):
    series.plot(label='Estimated population')
    decorate(xlabel='Year', 
             ylabel='Population estimate', 
             title='Narraguacus River',
             ylim=[0, 5000])
    
plot_population(pop_series)
Explanation: Here's what it looks like:
End of explanation
abs_diffs = pop_series.diff()
abs_diffs
Explanation: Modeling changes
To see how the population changes from year-to-year, I'll use diff to compute the absolute difference between each year and the next.
End of explanation
rel_diffs = abs_diffs / pop_series
rel_diffs
Explanation: We can compute relative differences by dividing by the original series elementwise.
End of explanation
rates = rel_diffs.dropna()
rates
Explanation: These relative differences are observed annual net growth rates.  So let's drop the 0 and save them.
End of explanation
np.random.choice(rates)
Explanation: A simple way to model this system is to draw a random value from this series of observed rates each year.  We can use the NumPy function choice to make a random choice from a series.
End of explanation
t_0 = 2015
p_0 = pop_series[t_0]
Explanation: Simulation
Now we can simulate the system by drawing random growth rates from the series of observed rates.
I'll start the simulation in 2015.
End of explanation
system = System(t_0=t_0,
                p_0=p_0,
                duration=10,
                rates=rates)
Explanation: I'll create a System object with variables t_0, p_0, rates, and duration=10 years. 
The series of observed rates is one big parameter of the model.
End of explanation
# Solution goes here
Explanation: Write an update functon that takes as parameters pop, t, and system.
It should choose a random growth rate, compute the change in population, and return the new population.
End of explanation
update_func1(p_0, t_0, system)
Explanation: Test your update function and run it a few times
End of explanation
def run_simulation(system, update_func):
    Simulate a queueing system.
    
    system: System object
    update_func: function object
    
    t_0 = system.t_0
    t_end = t_0 + system.duration
    
    results = TimeSeries()
    results[t_0] = system.p_0
    
    for t in linrange(t_0, t_end):
        results[t+1] = update_func(results[t], t, system)
    return results
Explanation: Here's a version of run_simulation that stores the results in a TimeSeries and returns it.
End of explanation
# Solution goes here
Explanation: Use run_simulation to run generate a prediction for the next 10 years.
The plot your prediction along with the original data.  Your prediction should pick up where the data leave off.
End of explanation
def plot_many_simulations(system, update_func, iters):
    Runs simulations and plots the results.
    
    system: System object
    update_func: function object
    iters: number of simulations to run
    
    for i in range(iters):
        results = run_simulation(system, update_func)
        results.plot(color='gray', label='_nolegend', 
                     linewidth=1, alpha=0.3)
Explanation: To get a sense of how much the results vary, we can run the model several times and plot all of the results.
End of explanation
# Solution goes here
Explanation: The plot option alpha=0.1 makes the lines semi-transparent, so they are darker where they overlap.
Run plot_many_simulations with your update function and iters=30.  Also plot the original data.
End of explanation
def run_many_simulations(system, update_func, iters):
    Runs simulations and report final populations.
    
    system: System object
    update_func: function object
    iters: number of simulations to run
    
    returns: series of final populations
    
    # FILL THIS IN
# Solution goes here
Explanation: The results are highly variable: according to this model, the population might continue to decline over the next 10 years, or it might recover and grow rapidly!
It's hard to say how seriously we should take this model.  There are many factors that influence salmon populations that are not included in the model.  For example, if the population starts to grow quickly, it might be limited by resource limits, predators, or fishing.  If the population starts to fall, humans might restrict fishing and stock the river with farmed fish.
So these results should probably not be considered useful predictions.  However, there might be something useful we can do, which is to estimate the probability that the population will increase or decrease in the next 10 years.  
Distribution of net changes
To describe the distribution of net changes, write a function called run_many_simulations that runs many simulations, saves the final populations in a ModSimSeries, and returns the ModSimSeries.
End of explanation
run_many_simulations(system, update_func1, 5)
Explanation: Test your function by running it with iters=5.
End of explanation
last_pops = run_many_simulations(system, update_func1, 1000)
last_pops.describe()
Explanation: Now we can run 1000 simulations and describe the distribution of the results.
End of explanation
net_changes = last_pops - p_0
net_changes.describe()
Explanation: If we substract off the initial population, we get the distribution of changes.
End of explanation
np.sum(net_changes > 0)
Explanation: The median is negative, which indicates that the population decreases more often than it increases.
We can be more specific by counting the number of runs where net_changes is positive.
End of explanation
np.mean(net_changes > 0)
Explanation: Or we can use mean to compute the fraction of runs where net_changes is positive.
End of explanation
np.mean(net_changes < 0)
Explanation: And here's the fraction where it's negative.
End of explanation
weights = linspace(0, 1, len(rates))
weights /= sum(weights)
weights
Explanation: So, based on observed past changes, this model predicts that the population is more likely to decrease than increase over the next 10 years, by about 2:1.
A refined model
There are a few ways we could improve the model.
It looks like there might be cyclic behavior in the past data, with a period of 4-5 years.  We could extend the model to include this effect.
Older data might not be as relevant for prediction as newer data, so we could give more weight to newer data.
The second option is easier to implement, so let's try it.
I'll use linspace to create an array of "weights" for the observed rates.  The probability that I choose each rate will be proportional to these weights.
The weights have to add up to 1, so I divide through by the total.
End of explanation
system.weights = weights
Explanation: I'll add the weights to the System object, since they are parameters of the model.
End of explanation
np.random.choice(system.rates, p=system.weights)
Explanation: We can pass these weights as a parameter to np.random.choice (see the documentation)
End of explanation
# Solution goes here
Explanation: Write an update function that takes the weights into account.
End of explanation
# Solution goes here
Explanation: Use plot_many_simulations to plot the results.
End of explanation
# Solution goes here
Explanation: Use run_many_simulations to collect the results and describe to summarize the distribution of net changes.
End of explanation
# Solution goes here
Explanation: Does the refined model have much effect on the probability of population decline?
End of explanation
<END_TASK>